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ABSTRACT 

The  immediate  attention  of  the  control  systems  engineer 
is  directed  to  the  dynamic  behavior  of  the  system  under 
study.   It  is  important  to  study  the  effects  on  overall 
system  performance  of  varying  one  or  more  parameters  (mass, 
inertia,  gain,  resistance,  etc.).   It  is  equally  important 
to  determine  whether  a  desired  dynamic  behavior  can  be 
achieved  with  any  set  of  values  for  the  parameters — if  not, 
redesign  is  indicated. 

In  this  thesis  a  control  systems  analysis  package  is 
developed  using  parameter  plane  methods.   It  is  an  inter- 
active, user-friendly  computer  aid.   Given  a  characteristic 
equation  containing  two  variable  parameters,  the  output  of 
the  analysis  may  be  either   tabular  or  graphical,  with  plots 
of  any  of  the  following  types: 

1.  Constant  damping  curves  as  a  function  of  frequency, 

2.  Constant  frequency  curves  as  a  function  of  damping, 

3.  Constant  sigma  lines  (real  root  lines), 

4.  Constant  zeta-omega  (damping-frequency)  curves. 
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I.   INTRODUCTION 

The  analysis  and  synthesis  of  linear  feedback  control 
systems,  or  the  compensation  of  same,  can  be  realized  by 
three  general  methods.   The  first  of  these  can  be  called  the 
integral  method.   Given  a  control  system,  described  by  a  set 
of  differential  equations,  One  selects  a  cost  function  to  be 
minimized  with  respect  to  certain  variable  system  parameters. 
The  major  drawback  with  this  method  is  the  difficulty  of 
varying  more  than  one  parameter  at  a  time.   The  second 
method  is  the  Bode  frequency  response  technique  whereby  the 
system's  open  loop  transfer  function  is  manipulated  to  obtain 
the  desired  system  response.   This  method  also  has  its 
inherent  weaknesses:   difficulty  of  application  to  non-unity 
feedback  control  systems,  difficulty  in  interpreting  the 
closed  loop  transient  response  in  terms  of  the  open  loop 
frequency  response,  and  difficulty  of  varying  more  than  one 
parameter.   Third  are  the  algebraic  methods.   Within  this 
category  can  be  included  the  familiar  root  locus  method. 
Here,  a  graphical  technique  is  provided  by  which  the  set  of 
all  points  which  could  potentially  be  made  roots  are  plotted 
in  the  S-plane.   The  root  locus  method  is  a  valuable  and 
powerful  tool  when  only  one  parameter  is  varied;  results 
are  less  satisfactory  for  two  parameters  and  of  little  use 
when  three  or  more  parameters  are  involved. 
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Methods  for  studying  the  parameter-root  relationship 
when  two  or  more  parameters  are  variable  are  clearly  of 
considerable  value.   For  a  linear  system,  the  set  of 
differential  equations  that  describe  that  system  can  be 
transformed  into  algebraic  equations  and  manipulated  to 
provide  a  characteristic  polynomial.   Since  the  coefficients 
of  the  characteristic  polynomial  are  deterimined  by  the 
system  parameters,  it  follows  that  some  relationship  exists 
between  the  value  of  any  parameter  and  the  value  of  the 
characteristic  roots.   In  reference  (1),  Mitrovic  developed 
an  algebraic/graphical  method  for  obtaining  the  roots  of  a 
polynomial  in  terms  of  two  variable  parameters.   In 
references  (2),  (3),  and  (4),  Choe ,  Hyon ,  and  Nutting, 
respectively  developed  and  extended  the  Mitrovic  method  to 
the  compensation  of  linear  continuous  feedback  control 
systems.   The  disadvantage  of  the  Mitrovic  method  is  that  the 
variable  parameters  may  appear  in  no  more  than  two  coeffi- 
cients of  the  characteristic  equation,  which  limits  the 
flexibility  of  the  technique.   In  reference  (5),  Siljak 
introduced  a  method  for  obtaining  the  roots  of  a  polynomial 
in  terms  of  two  variable  parameters  that  may  appear  in  any 
and  all  of  the  coefficients  of  the  polynomial.   Later,  Thaler 
and  Towill  [Ref.  6]  extended  this  method  to  the  compensation 
of  linear  continuous  feedback  control  systems.   It  is  from 
the  latter  work  that  the  ensuing  parameter  plane  equations 


were  developed.   General  methods  of  compensation  will  be 
presented,  and  an  attempt  will  be  made  to  relate  the  root 
locus  and  the  parameter  plane  methods  as  a  set  of  comple- 
mentary techniques  which,  when  applied  in  tandem,  represent 
the  most  satisfactory  tool  to  date  for  designing  linear 
feedback  control  systems. 

The  parameter  plane  method,  which  works  well  for  two 
variable  parameters  and  which  may  be  extended  to  three  or 
more  parameters,  is  purely  algebraic,  and  the  resulting  plots 
are  valuable  aids  to  analysis.   The  term  parameter  plane 
comes  from  the  plot  for  two  parameters  —  in  a  rectangular 
coordinate  space  one  parameter  will  define  the  abscissa  while 
the  second  parameter  defines  the  ordinate  (the  S-plane  is 
inconvenient  for  presenting  the  desired  results).   Three 
parameters  define  a  3-dimensional  parameter  space,  etc.   For 
design  problems  it  is  convenient  to  think  of  the  algebraic 
calculations  as  a  mapping  procedure.   By  choosing  a  point  on 
the  S-plane,  the  characteristic  polynomial  acts  as  a  mapping 
function  whereby  the  point  may  be  "mapped"  onto  the  alpha- 
beta  plane  (alpha  and  beta  are  the  two  variable  parameters 
to  be  used   throughout  the  remainder  of  this  text).   The 
relationship  between  being  able  to  place  the  roots  of  a 
polynomial  at  specific  locations  in  the  S-plane  and  the 
compensation  of  linear  feedback  control  systems  is  as 
follows.   A  feedback  control  system,  including  any  added 


compensators  which  may  contain  variables,  can  be  reduced  to  a 
ratio  of  two  polynomials  (the  closed  loop  transfer  function). 
A  specified  system  response  in  terms  of  overshoot,  bandwidth, 
settling  time,  steady-state  accuracy,  etc.,  can  theoretically 
be  obtained  by  placing  a  pair  of  complex  conjugate  roots  of 
the  characteristic  equation  at  a  specific  location  in  the  S- 
plane,  while  ensuring  that  the  real  part  of  this  complex  root 
pair  (the  dominant  roots)  is  smaller  in  magnitude  than  the 
real  parts  of  the  remaining  roots  of  the  characteristic 
equation.   The  problem  of  compensation,  thus  of  feedback 
control  system  design,  reduces  to  one  of  placing  the  dominant 
roots  of  the  characteristic  equation  at  the  desired  location. 
The  ability  of  the  parameter  plane  method  to  achieve  this 
goal  will  become  obvious. 


II.   DERIVATION  OF  PARAMETER  PLANE  EQUATIONS 

A  linear  feedback  control  system's  characteristic  equation 
can  be  expressed  as  a  polynomial  of  the  following  form: 


m     k 
f(s)  =  Z      a,S   =  0,  where  (2-1) 

k=0   K 


a,  (k=0 ,1 , . . . ,m)  are  real  coefficients 
S  =  -a+jw  =  -|cu+  3<*i   Jl-£2 
w  is  the  undamped  natural  frequency  and 
?^  is  the  relative  damping  coefficient 

In  reference  (5)  it  is  noted  that  S   may  be  represented  by 

the  following: 


Sk  =  cjk(Tk(-|)  +  j  /  1-|2  U  (-O)  (2-2) 


where 

Tk(-C)  =  (-DkTk(0  and  Ufe(  -|)-(  -1  )k+1Uk(  |)  . 

T  (^)  and  U  (£)  are  Chebishev  functions  of  the  first  and 
second  kind  respectively.   Values  of  zeta  and  omega  will  be 

considered  such  that  0<£<1  and  0<cu<oo.   Values  of  T.  and  U, 

— ^  —       —  —  k       k 

are  tabulated  in  -various  appendixes.  More  important  to  digital 
computer  analysis,  they  can  be  obtained  from  the  following 
recursive  relations: 
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Tk+ia)  -  2Tk(*}  +  Vi(n  -  ° 


W*>  -  aV*>  +  uk-i(^  =  ° 


(2-3) 


Here,  TQ(|)=1,  1^(1  )-£,  UQU)=0,  1^(0  =  1.   Substituting 
equation  (2-2)  into  (2-1)  and  setting  the  real  and  imaginary 
parts  to  zero  independently,  one  obtains: 
m 


2   a  cukT  (-0  =  0 
k=0   K    K 


m     k 

Z   a,cAj,(-£)=0 

k=0   K    K 


(2-4) 


Employing  equations  (2-3),  one  obtains  from  equation  (2-4) 


m  lr     lr 

2   (-DVcAj   -U)  -  0 
k=0       K    K_1 


(2-5) 


m  lr     lr 

I   (-1)  a.«\(5)  =  0 
k=0       K    K 


Now  consider  the  coefficients  a.  of  the  characteristic 

k 

equation  (2-1)  as  linear  functions  of  the  variable  system 
parameters,  a   and  /3,  as  follows: 


ak  ■  V  +  V3  +  dk  (2-6) 
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Using   this    relation    for    a,  ,    equations  (2-5 )    become: 


aB     +   3C      +   D      =    0 


(2-7) 


°'B2    +5C2    +    D2    =    ° 


where 


k   ,       k   TT  „  "f      ,    _k   ,       k 


B_    =      S      (-1)      b,a)      U,  B      =      S      (-1)      b,a>     U. 

1         k=o  K  K_1  z         k=0  K  K 


C      =      S      (-1)K   ck«*  U  C      =      2      (-1)K   c   «,*  U 

1  k=0  K  K     1  Z  k=Q  K  K 


(2-8) 


m  ,  ,  m  .  , 

D      =      2      (-1)*   d  o>*   U  D      =      2       (-1)K   d  u*   (J 

1        k=0  k  k-i  ^        k=0  K  R 

Since   equations  (2-7 )  are  linear    in    the    two    unknowns      alpha    and 
beta,    Cramer's    rule   may   be    applied   to   obtain: 


C1D'>-C2D1  B2D1"B1D2 

BlVB2Cl  B1C2-B2C1 
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Equations (2-9)  are  now  functions  of  zeta  and  omega.   Hence, 
by  fixing  either  zeta  or  omega  and  varying  the  remaining 
parameter,  the  constant  omega  or  constant  zeta  S-plane 
contours  respectively  can  be  mapped  into  the  real  domain 
of  the  alpha-beta  or  parameter  plane. 

In  reference  (5)  the  following  relationships  are  noted: 


SK  =  Pk  +  ju/  1-^      Qk 


Pk+1  +  2UPk  +  U2pk-1  "  ° 


Qk+1  +  2uQk  +  u2Qk_1  -  0  (2-10) 


Pk  =  -,^.Qk  -  o,2Qk_1 


Po  =  %   =  °.  Pl  =  -«u>  Ql  =  l 


P,  and  Q,  are  related  to  the  Chebishev  functions  by: 


Pk  =  ^V-0  =  (-Dka'kTk(0 


Qk  =  u^^c-n  =  (-Dk+1a)k-1uk(u 


(2-11) 


12 


By  employing  equations  (2-10)  and  (2-11),  one  obtains 
(proceeding  as  before); 


m 

z   akQk-l  =  ° 
k=0   K  K  i 


m 

2   a  Q   =  0 
k=0   K  K 


(2-12) 


Combining  equations  (2-6)  and  (2-12)  with  Cramer's  rule,  one 
again  arrives  at  equations  (2-9)  where  the  following  expressions 
now  apply: 


m 

Bl  =   S   bk  Qk-1  =  ° 
1    k=o   K   K  x 


m 


B2  =   Z  \   Qk  =  ° 
A         k=0   k   K 


Cl  =  Z  ck  Qk-l  =  ° 
1    k=0  K      K  L 


m 

C2  =   2-   Ck  Qk  =  ° 
Z    k=0   K   K 


(2-13) 


D. 


m 
=   I 
k=0 


dk  Qk-1  "  ° 


m 

D2  =   2   dk  Qk  =  ° 
Z         k=0   R   k 


Equations  (2-9)  and  (2-13)  are  useful  for  mapping  constant 
zeta-omega  curves  from  the  S-plane  into  the  parameter  plane. 
As  will  be  demonstrated  later,  these  curves  play  an  important 
role  in  dominance  considerations. 

If  the  complex  variable  S  is  substituted  in  equation  (2-1) 
by  letting  S  =— a,  where  sigma  corresponds  to  values  of  S  along 
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the  real  axis,  then  according  to  equation  (2-6)  the 
characteristic  equation  (2-1)  becomes: 


a   Z      (-D  K°      +     P    2   (-1)  c  a   +  S   (-1) Kd  °K  =  o 
k=0       K       k=0       K     k=0       K 


(2-14) 

The  above  expression  represents  a  straight  line  in  the  alpha- 
beta  plane  for  a  given  value  of  sigma.   Hence  a  point  on  the 
real  axis  in  the  S-plane  maps  into  a  straight  line  in  the 
alpha-beta  plane.   In  addition,  for  given  values  of  alpha, 
beta,  and  sigma  which  satisfy  equation  (2-14),  the 
characteristic  equation  (2-1)  must  have  a  real  root  at  minus 
sigma.   On   the  constant  zeta  and  omega  curves  previously 
defined,   for  certain  values  of  alpha  and  beta  (say,  for 
values  obtained  from  equations  (2-9)  for  given  values  of 

zeta  and  omega)  the  characteristic  equation  will  have  a  pair 


2 
of  complex  conjugate  roots  at  S  =  -§a;  +  joi/  1-C 

The  significance  of  the  above  discussion  is  that  by 

applying  equations  (2-9)  and  (2-14)  one  can,  for  a  specified 

value  of  zeta,  omega,  and  sigma,  compute  the  value  of  alpha 

and  beta  such  that  the  characteristic  equation  will  have  a 


2 
pair  of  roots  at  S  =  -£-,<*>-,  +  jw-V  l-£,    .   The  m-2 

remaining  roots  of  the  characteristic  equation  can  then  be 

calculated  by  dividing  out  the  two  known  or  specified 

roots.   This  method,  where  zeta,  omega  and  sigma,  or  simply 
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zeta  and  omega  are  specified,  and  where  the  computations  for 
alpha  and  beta  are  done  algebraically,  will  be  referred  to  as 
the  algebraic  parameter  plane  solution. 

To  solve  the  problem  in  general  for  all  values  of  zeta, 
omega,  and  sigma,  it  becomes  necessary  to  plot  a  family  of 
parameter  plane  curves  for  various  values  of  zeta,  omega, 
sigma,  and  if  desired,  zeta-omega.   On  the  resulting 
parameter  plane  plot  one  can,  by  choosing  an  operating  point, 
graphically  read  from  the  curves  the  values  of  alpha  and  beta 
and,  hence,  the  values  of  the  m  roots  of  the  corresponding 
m   order  characteristic  equation.   This  latter  method  will 
be  referred  to  as  the  graphical  parameter  plane  solution. 

The  algebraic  solution  has  the  advantage  that  the  labor 
of  plotting  the  curves  can  be  avoided,  but  the  disadvantage 
remains  that  without  the  curves  it  is  difficult  to  pick  the 
optimum  values  of  zeta  and  omega  so  as  to  ensure  dominance 
while  still  meeting  the  system  specifications.   The  graphical 
solution  has  the  advantage  that  one  has  a  "picture"  of  the 
way  the  characteristic  roots  move  about  in  the  S-plane  as 
alpha  and  beta  are  varied.   This  enables  one  to  choose  the 
values  of  alpha  and  beta  corresponding  to  the  best  values  of 
zeta,  omega,  sigma,  and  zeta-omega  for  all  roots  of  the 
characteristic  equation.   This  feature  of  the  parameter 
plane  points  out  a  strong  justification  for  attempting  to 
obtain  the  parameter  plane  curves.   And  with  the  employment 
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of  a  digital  computer  and  an  appropriate  algorithm  to  realize 
the  parameter  plane  curves,  the  advantage  of  the  algebraic 
method  becomes  muted. 

Recursion  methods  (equation  (2-3))  are  by  no  means  the 
only  methods  of  producing  algebraic  and  graphical  parameter 
plane  data.   Thaler  and  Karmarkar   [Ref.  7]  describe  a  matrix 
solution  to  the  parameter  plane  problem.   Essentially,  a 
matrix  of  coefficients  may  be  manipulated  to  obtain  the 
following  general  form: 


b        c.       -w 


b. 


0 


-2£u>      -u) 


-1 


-2$w 


0 


0 


•a> 


-1  -2$u> 


b         c         0  0  0 

m        m 


0 


0 


0 


-«*> 


m-2 


-2*>      d      .-a* 
m-1 


-1  d    -2tw        e 

m      ? 


m 


m- 


R 


m-2 


m-2 
R. 


a8 


=  0 
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where  b,  ,  c  ,  d  ,  4 ,  and  gwu  are  as  described  before,  and: 


m)  are  the  coefficients  associated  with  the 


ek  (k=l,  . . 

non-linear  alpha-beta  product  terms 


m-^ 


R1  is  the  sum  of  the  m-2  roots  of  the  polynomial  character- 
istic equation  taken  one  at  a  time 


m-: 


R   0  is  the  sum  of  the  m-2  roots  taken  m-2  at  a  time 
m-2 


Further,  by  the  application  of  appropriate  row  operations, 
this  matrix  may  be  reduced  to  the  following  row-echelon  form: 


10        0 


0        10 


0        0        1 


0         0         0 


0      kn      k12 


0         k21         k22 


0         k31         k32 


1         k   ,         k   _ 
ml  m2 


m-2 
R 


m-2 


aB 


=    0 


(2-15) 
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For  the  case  when  all  coefficients  of  the  characteristic 
equation  are  linear,  i.e.  e,  (k=l ,  .  .  .  ,m)  =  0,  then  K 
(k=l, . . . ,m)  =  0,  and 


iVll 

6  =  -KOH 
21 

as  obtained  from  the  first  two  rows  of  the  matrix  equation. 

One  should  note  that  in  arriying  at  equations  (2-15)  , 

2 
approximately  m   row  operations  are  required  for  the  row- 
echelon  matrix  formulation  for  each  point  of  the  parameter 
plane  curves  (e.g.,  each  time  either  zeta  or  sigma  are 
varied).   Compare  this  with  the  approximate  m  calculations 
required  to  obtain  the  recursion  equations  of  the  previous 
chapter,  and  the  matrix  method  becomes  relatively  inefficient 
for  larger  order  systems. 

One  should  not,  however,  discard  the  matrix  approach 
entirely.   For  small  order  characteristic  equations,  this 
technique  compares  favorably  with  the  recursion  method. 
And  when  the  variable  parameters  are  non-linear--when  one 
must  deal  with  alpha-beta  product  terms—the  matrix  approach 
affords  a  more  direct  method  of  obtaining  the  alpha-beta 
pairs.    Whether  the  recursion  or  matrix  method  is  utilized 


From  equations  (2-15),  one  obtains  the  two  quadratic  forms 

K22a'-(K21  K12-Kll  K22-1)a+Kll=° 

K126  +(K21  K12-Kll  K22+1)6+K21=° 
from  which  alpha  and  beta  are  easily  derived. 
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should  depend  on  the  inclusion  of  alpha-beta  product  terms; 
ultimately,  it  is  a  matter  of  personal  preference. 
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III.   APPLICATION  OF  THE  PARAMETER  PLANE  METHOD 

A.   ALGEBRAIC  SOLUTION 

In  this  section  it  will  be  assumed  that  the  system 
performance  specifications  can  be  met  by  placing  a  pair  of 
complex  conjugate  roots  at  a  specific  location  (i.e.,  by 
choosing  appropriate  values  for  zeta  and  omega).   If  after 
computation  of  the  necessary  values  of  alpha  and  beta  to 
locate  the  roots  as  desired  it  is  found  that  these  specified 
roots  are  not  dominant ,  then  either  a  different  value  of  zeta 
and/or  omega  must  be  used  (possibly  at  the  sacrifice  of  some 
performance  measure),  or  a  different  method  of  compensation 
will  have  to  be  attempted.   In  a  later  section  a  method  will 
be  addressed  whereby  the  dominancy  requirement  may  be. 
achieved. 

1.   Feedback  Compensation 

For  a  unity  feedback  control  system,  let 


q  _   K    _  K (3-1) 

e(S)    Sm+e   1Sm-l+...+eTSL 
m-1  L 

where  K  is  the  forward  path  gain  (a  variable)  and  e(S)  is  a 
polynomial  in  S  representing  the  poles  of  the  open  loop 
transfer  function  of  the  uncompensated  system.   In  equation 


20 


(3-1),  L  corresponds  to  the  system  type—for  a  type  0 
system,  L=0,  for  type  1,  L=l,  etc.   The  system's  error 
coefficient  is  defined  as: 


K   =  lim   SLG  (3-2) 

e    s^0      cc 


where  G    is  t'he  open  loop  transfer  function  of  the 

CC  f  f 

compensated  system. 

a.   Tachometer  Plus  Acceleration  Feedback 

In  order  to  achieve  the  system  performance 
specifications,  a  feedback  compensator  must  be  introduced. 
Let 

H  =  K  S  +  K  S2 
t      a 

The  resulting  compensated  system's   characteristic  equation 
becomes: 

e(S)  +  K(K.S+K  S2)  =  0  (3-3) 

l    a 

and  by  expanding  e(S),  equation  (3-3)  becomes: 

Sm+e   nSra_1+. . .+(e  +KK  )S2+(  e., +KK  )S+e  +K  =  0      (3-4) 
m-1  2    alto 

where  L  is  zero  for  a  type  0  system  (the  most  general  case). 
The  following  results  also  apply  to  a  type  1  system  if  e   is 
set  to  zero,  and,  similarly,  for  a  type  2  system  if  both 
eQ  and  e,  are  set  to  zero,  etc.   Combining  equations  (3-2), 
(3-3),  and  (3-4)  the  error  coefficient  becomes: 

21 


K   =  lim  5-1 =  JL-  (3-5) 

e    S+0   e(S)+K(K  S+K  S  )     o 

t    a 

or  for  a  type  1  uncompensated  system: 

K   =    \v„  (3-6) 

e    e  ,+KK 

or  for  a  type  2  uncompensated  system: 

Ke  ■  n;  .  (3-7) 


Note  that  if  the  uncompensated  system  is  type  2,  the 
compensated  system  would  be  type  1  if  tachometer  feedback 
or  tachometer  plus  acceleration  feedback  is  used. 

In  the  compensated  system's  characteristic 
equation  (3-4)  let  alpha  =  KK   and  beta  =  KK  .   Equation 
(3-4)  then  becomes: 

Sm  +  e  ,  S"1"1  +.  .  .  +  (e0  +  ct)S2+(e _  +  8)S+e.  +K  =  0 
m-i  a  1      o 

Recalling  equation  (2-6)  where  in  general  the  coefficients 
of  the  characteristic  equation  are  of  the  form: 


ak  =  V+ckB+dk 


and  letting  m=k,  then  from  equations  (2-8)  one  obtains: 
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B   =  (-1)-^2U1  =  "2  B2  =  "2U2 

C,  =  -wU   =0  C0.=  -«U_  =  -a)  (3-8) 

1       o  2       1 

m       k    k  m>      k    k 

D   =   E   (-D  VX-1  °2  =  Z       (-1}  V  Uk 

1     k=Q  K     K  1       Z     k=0  K     k 


since  U  =0  and  U  =1.   From  equations  (2-9)  one  derives: 
o        1 

m      k 

C1D2-C2D1     k=0  (_1)  dkUk"1      m       k    k-2 

1  J Z_J^_    _        K-U (-l)Kd  °»     U 

12   2  1  -co  k=0 


(3-9) 


k=0 


If  alpha  and  beta  are  linear  functions  of  K,  the  forward 
path  gain,  one  can  use  the  steady-state  error  specification 
to  define  K  in  terms  of  alpha  and/or  beta.   Since  zeta  and 
omega  were  assumed   to  be  specified,  then  from  equations 
(3-9)  one  may  solve  for  alpha  and  beta.   From  this,  K   and 

K^  are  readily  determined. 

t  J 

Example  3-1 
The  system  of  Figure  (3-1)  is  to  be  compensated  by  using 
tachometer  plus  acceleration  feedback.   The  system 
specifications  are  as  follows: 
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%£>. 


A 


S(S+1)(S+2) 


KQS  +  KtS 


<- 


-> 


Figure  3-1 


1.  Complex  roots  corresponding  to  £  =  0.7  and  «>=io 

2.  K  >6 

e— 


From   equation    (3-2) 


K      = 


K 


>  +  KK 


>    6 


From   which    K_>12+6KK    .       The    compensated    characteristic 
equation    is 


S3  +  S2(3  +  KK    )  +  S(2  +  KK4.)+K   =    0 
a  t 


(3-10) 
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Letting  alpha  =  KK   and  beta  =  KK  ,  equation  (3-10)  becomes: 

3,  L 


S3+S2(3+a;)  +  S(2  +  i3)+K  =  0  (3-10a) 


From  equations  (3-8): 


B   =  100  B2  -    140 

C   =  0  c2  =  -10 

D   =  -1100-K  D   =  -1120 


and  from  equations  (3-9): 


„    _    lO(-llOO-K)      a   _  140(-1100-K)+56000      ,  „    . 

^Tooo    '  p  =Tooo         ^-11; 


From  the  steady-state  accuracy  specifications,  then,  it  is 

necessary  that  K>12+6/3;  let  K=12+6/3.   From  equations  (3-11) 

it  is  found  that  0=62  3,  hence  K=3750.   Therefore,  a=48.5, 

and  since  a=KK   and  /3=KK,  : 
a         t 


48    F5 
K      =   |§~   =    0.0129 
a         3750 


Kt  =  §§o  "  °-1661 
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The  compensated  system's  characteristic  equation  becomes 


S3+51.5S2+625S+3750=0  (3-12) 


2 
Now  zeta=0.7  and  oo  =  10  corresponds  to  S  +14S+100=0. 

Dividing  equation  (3-12)  by  this  quadratic,  the  remainder 

is  S+37.5.   Since  zeta*omega  of  the  desired  roots  =  7<<37.5, 

the  complex  roots  are  dominant  and  the  problem  is  solved. 

b.   Tachometer  Feedback  Only 

Let  H  =  K  S.   The  characteristic  equation  of  the 

compensated  system  becomes: 


Sm+em_1Sm  1+.. .+e2S2+(e1+a)S+e0+6=0 


Proceeding  as  in  the  previous  example,  one  obtains 


B1    =  0  B2  =  -a) 


Cl  "  -1  C2  =  ° 


m        k      k  m        k    k 

D-  =   Z    (-D  d.   coKU,         D   -   z    (-l)\u\ 
1   k=0        k     K-i       d        k=0        k   k 
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and 


m       ,    .  ,  m       ,    , 

a  =  Z   (-1)VK  V  ,     3  =  Z   (-1)V\  , 
k=0        K      K  k=0        K    K_i 


(3-13) 


For  a  specified  value  of  zeta  and  omega,  alpha  and  beta  can 
be  obtained  from  equations  (3-13).   The  error  coefficient  is 
then  determined  directly  from  equations  (3-5),  (3-6),  or 
(3-7).   Thus  the  error  coefficient  is  fixed  for  a  given 
value  of  zeta  and  omega,  and  if  this  parameter  is  to  be  met, 
the  values  of  zeta  and  omega  may  require   adjustment.   One 
possible  approach  might  be  to  fix  zeta  at  some  value, 
whereby  from  the  given  K   and  equations  (3-5),  (3-6),  or 
(3-7)  alpha  could  be  computed.   Equations  (3-13)  could  then 
be  solved  for,  first,  omega  and  then  beta.   The  calculations 
would  prove  tedious,  however. 
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Example  3-2 
The  same  system  as  used  in  example  (3-1)  will  be  studied 
here,   this  time  with  tachometer  feedback  alone.   The  same 
system  performance  specifications  are  to  be  met,  namely, 
K  >_6 ,  zeta  =  0.7,  and  omega  =  10.   The  compensated  system's 
characteristic  equation  becomes: 


S3+3S2+(2+KK  )S+K  =  0  (3-14) 


Letting  a=KK   and  /3=  K  here,  equation  (3-14)  becomes 


S3+3S2(2+a)S+/3  =  0 


From  equation  (3-13)  it  is  found  that: 

a   =  -2+30(1. 4)-100(0. 96)  =  -56 

Since  alpha  is  negative,  it  is  seen  that  positive  tachometer 
feedback  is  required.   Further,  it  is  found  that  the 

remaining  root  (when  equation  (3-14)  is  divided  by 

2 
S  +14S+100)  is  positive;  the  system  is  unstable.   Hence 

the  desired  system  specifications  cannot  be  met  with 

tachometer  feedback  alone. 

c.   Acceleration  Feedback  Only 

2 
Let  H  =  K  S  .   The  characteristic  equation  of 
a 

the  compensated  system  becomes: 
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Sm+e   ,Sm  1+...+(e„+KK  )S2  +  e.,S+e  +K  =  0 
m-1  2    a      1    o 


Proceeding  as  before,  where  now  a=KK   and  fS=K 

a 


B   =  -2V    =»2 


B2  "  "  U2 


Cl  =  U-l  "  -1 


C2  "  ° 


m      ,    , 
D   =  2   (-l)Kd  k-X-l 

1     k=0  K     K  1 


m      ,    , 

D   =   Z   (-l)Kd  »*U 
*         k=0        K    K 


Solving  for  alpha  and  beta  yields: 


-D, 


a  = 


u2U, 


—  (-l)kd  "k~2U 
U2  k  ±;  ak     uk 


m      k    k       1         m  k    k 

0.   X   (-l)kd  - ku     -  i-  2   (-l)kdokU 
k=0         K    k  i     U2  k=Q         k    k 


(3-15) 


Calculations  for  alpha,  beta,  and  K   are  performed  in  the 
same  manner  as  with  the   preceding  tachometer  feedback 
example . 

Example  3-3 
The  same  system  of  examples  (3-1)  and  (3-2)  will  now  be 
compensated  using  acceleration  feedback  alone.   As  before, 
K  >6 ,  zeta=0.7,  and  omega=10.   Therefore  K  =~  and  the  error 

*3  t?   ju 
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coefficient  is  unaffected  by  the  acceleration  feedback.   Hence 
one  can  conveniently  choose  K=12  to  meet  the  specifications. 
The  compensated  system's  characteristic  equation  becomes: 


S3+(3+KK  )S2+2S+K  =  0  (3-16) 

a 


If  K  in  equation  (3-16)  is  set  equal  to  12  as  prescribed,  only 
one  parameter  remains  and  the  parameter  plane  equations 
produce  an  indeterminate  solution.   If  K  is  left  as  the 
variable  beta,  then  equation  (3-16)  becomes  (after  the  usual 
subst itut  ions ) : 


S3+(3+a)S2+2S+/3  =  0 


By  employing  equations  (3-15),  one  obtains  a=4  and  j3  =  -700. 
Since  beta  is  negative,  it  is  concluded  that  the  desired 
roots  (i.e.,  desired  values  of  zeta  and  omega)  cannot  be 
realized  using  acceleration  feedback  alone,  and  of  course 
neither  can  the  desired  error  specification   be  obtained. 
One  would  therefore  choose  an  alternate  method  of 
compensation . 

If  one  chooses  to   use  feedback  compensation  then 
perhaps  tachometer  plus  acceleration  feedback  might  be 
attempted  first  using  equations  ( 3-9 )  and  the  appropriate 
steady-state  error  specification.   If  the  specif ications  cannot 
be  met  in  this  manner,  then  it  follows  that  neither 
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tachometer  nor  acceleration  feedback  alone  will  suffice.   In 

this  case  either  the  system's  specifications  must  be  cased  or 

another  type  of  compensation  must  be  utilized.   If  it  is 

found  that  the  specifications  are  achievable  with  the  combined 

tachometer  and  acceleration  feedback,  then,  if  desired, 

equations  (3-13)  and  (3-15)  can  be  employed  to  investigate 

the  feasibility  of  tachometer  or  acceleration  feedback  alone, 

respect  i ve  ly . 

d.   Case  For  Which  Feedback  Is  Not  Available  Near 
The  Forward  Path  Amplifier 


+/^~N 


H+f-M 


+ 


1 


S(S+1)(S+-2) 


KQS  +  KtS 


Figure  3-2 

Figure  3-2  shows  a  system  similar  to  that  used 
in  example  (3-1)  except  that  now  the  feedback  is  inserted  at 
the  output  terminals  of  the  amplifier  represented  by  gain  K. 
This  illustrates  a  system  for  which  it  may  not  be  possible 
or  practical  to  access  the  input  terminals  of  the  error 
detector.   This  problem  will  be  solved  by  means  of  an  example 
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Example  3-4 
As  before,  the  same  system  specifications  are  to  be  met, 
i.e.,  K  >6,    zeta=0.7,  and  omega=10.   The  characteristic 
equation  becomes: 


S3+(3+K  )S2+(2+K  )S+K  =  0 
a        t 


Letting  a=K   and  0=K  : 

Si  X 


S3  +  (3+a;)S2  +  (2  +  /3)S+K  =  0  (3-17) 


Comparison  of  equations  (3-17)  and  (3-10a)  show  that  they 
are  identical,  that  is,  the  solution  obtained  for  alpha  and 
beta  in  example  (3-1)  applies.   There,  alpha  and  beta  were 
found  to  be  48.5  and  623,  respectively,  while  K=3750.   For 
the  present  example  no  further  computations  are  necessary 
to  find  K   and  K  ,  since  they  are  now  the  parameters  alpha 
and  beta.   This  points  out  an  important  advantage  of  the 
parameter  plane  method,  namely,  that  the  solutions  depend  only 
on  the  characteristic  equation  and  not  on  the  system  from 
which  the  characteristic  equation  was  formed.   This  principle 
can  similarly  be  applied  to  control  problems  involving 
tachometer  or  acceleration  feedback  alone. 
2 .   Cascade  Compensation 

For  a  unity  feedback  control  system  let  G  have  the 
form  of  equation  (3-1): 
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G  =  -T-T-  =  —, r  (3-1) 

e(s)    sm+e    sm-l+!>>+   SL 

m-1  L 

where  again  K  is  the  forward  path  gain  (a  variable)  and 
e(s)  is  a  polynomial  in  S  representing  the  poles  of  the 
open  loop  transfer  function  of  the  uncompensated  system. 
The  letter  L  again  indicates  the  system  type.   If  in  order 
to  satisfy  the  system's  requirements  a  cascade  compensator 
Gp  is  required,  then  let: 


P(S+Z) 
UC    Z(S+P) 


With  a  d.c.  gain  of  unity,  placement  of  this  compensator  in 
the  forward  path  will  not  affect  steady-state  accuracy. 
With  Gc  as  indicated  here,  the  values  of  P  and  Z  are 
computed  to  obtain  the  desired  system  response.   If  P  is 

less  than  Z,  a  lag  network  is  required  and  the  factor  of 

p 

■y   of  the  compensator  is  inherently  present  due  to  the 

physical  nature  of  the  compensator  (usually  an  R-C  network) 

In  this  case  all  forward  path  amplifier  gains  can  remain 

unaltered  to  meet  the  specified  accuracy  demands.   If, 

however,  the  computed  value  of  P  is  greater  than  Z,  a 

lead  network  is  required  and  the  compensated  system's 

p 
forward  path  gain  must  be  raised  by  the  factor  of  y   to 

meet  the  accuracy  specifications.   As  the  physical  nature 

p 
of  the  lead  network  is  such  that  the  factor  y    is  not 
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inherently  present,  this  factor  must  be  provided  either  by 
adding  an  amplifier  in  cascade  with  the  lead  network  or  by 
raising  the  gain  K  of  the  existing  amplifier  as  required  to 
achieve  steady-state  accuracy. 

Continuing,  the  compensated  system's  forward  path 
transfer  function  is: 

p 
G     G  G  =   K     P    S+Z    7(S+y)     K 


CC   "C     e(s)    Z    S+P     S+P     e(s) 

Applying  the  definition  of  the  error  coefficient  one 
obtains: 


is-  i    ■  cX'       r   K  Y     -.  K 


and  again  assuming  a  type  0  system  where  L=0,  the  compensated 
system's  characteristic  equation  becomes: 


e(s)(S+P)+KY(S+-)  =  0 


or  after  expansion: 


Sm+1  +  (P+em  ,)SmHPem   ,+em  2)Sm-1  + 

m-1        m-i   m-z 


+(Pen+en )S2+(KY+Pe  +e  )S+P(e  +K)  =  0  (3-18) 

2.       1  1   o       o 
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Letting  ot  =p    and  3  =  y,    equation    (3-18)    becomes: 


Sra+1  +  (a  +  em    .)Sm+(em     a+e      ^s"1"1*... 
m-1  m-i        m-6 


+(e2a+e1)S2+(K3+eia+e0)S+a(e0+K)  =  0  (3-19) 


Comparison  of  equation  (3-19)  with  the  general  form  of  the 
characteristic  equation  as  specified  in  equations  (2-1)  and 
(2-6),  it  is  apparent  that  K=m+1,  b^e^+K,  co=do=0,  b  =  e ,  , 

c  =K,  d,=eQ,  b2=e2 ,  c2=0>  d2=el'  etc" 

It  is  important  to  note  that  the  parameter  plane 
variable  beta  represents  the  pole-to-zero  ratio  of  the 
cascade  compensator.   The  S-plane  can  be  divided  into  regions 
where  lag  compensation  or  lead  compensation  is  needed.   By 
mapping  of  variables  in  the  above  manner,  the  parameter 
plane  can  effectively  be  divided  into  corresponding  regions 
above  and  below  the  straight  line  3=1.   Then,  for  values  of 
beta  less  than  one  a  lag  network  is  required  and  for  beta 
greater  than  one  a  lead  network  is  needed.   In  addition,  if 
beta  is  less  than  0.1  or  greater  than  10,  a  multiple  lag  or 
multiple  lead  network,  respectively,  is  required. 

Based  on  equations  (3-19)  and  (2-8)  it  is  found  that: 


Bx  =  -(e0+K)  +  a32e2+...+(-l)k"2^k-2Uk_3+(-l)k-1o)k  1Uk_2 


Cl  =  ° 
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Dn  =  o,2e,  +  .  ..+(-l)k"2e   Q-k_2U1  _  +  (-i)k-1e   -k"1U1  _ 
1       1  m-2      k-3  m-1      k-2 


+  (-Dku,kuk_1 


B2  =  -.e1^2e2U2+...+(-l)k-2a)k-2Uk_2  +  (-l)k-1,k-1Uk_i 


C2  =  -o»K 


D2  ,  -^o+.2e1U2+...  +  (-l)k-2em_2^k-2Uk_2 


+  (-l)k-1e   ,«k"1U1   +(-l)k"'kU1 
m-1      k-1  k 


(3-20) 


and  from  equations  (2-9)  it  is  found  that: 


„    "Dl      .        B2D1-B1D2        •  f,  „, 

a  =   p —   ,  P  =   p— 5 (3-21) 

Bl  C2B1 


For  a  type  1  system,  e   in  equations  (3-20)  is  set  equal  to 

zero,  for  a  type  2  system  e  =e  =0,  etc.   On  the  basis  of 

equations  (3-20)  and  (3-21)  a  cascade  compensator  can  be 

designed. 

Example  3-5 

For  the  system  of  Figure  (3-3)  it  is  desired  to  design  a 

cascade  compensator  which  places  a  pair  of  roots  at  zeta=0.5 

and  omega=l.   The  error  coefficient  K   should  be  50. 
&  e 
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+  h-^ 


(S+P/  ) 

S  +  P 


;S+1)(S+2)(S  +  b 


Figure    3-3 

1/ 

It  is  apparent  from  Figure  (3-3)  that  K=—j=50,  or 
K=500.   The  characteristic  equation  is: 


S4+(8+P)S3+(l7+8P)S2+( 10+17P+KV)S+P( 10+K)=0 


and  by  the  substitutions  of  *  =P  and  3=  Y: 


S4+(8+aOS3+(17  +  8LOS2+(10  +  17^  +  500£)S  +  510U'   =    0 


Applying  equations  (3-20)  and  (3-21)  one  obtains 


B. 


D,  = 


a      = 


-503 

0 

9 

0.0179 


=  P 


B2  = 
C2  = 

D2  = 

f      = 


-9 

-500 
6 
0.0117 


=  y 
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p 
and  since  7 '  =  %,    Z=1.529.   This  is  a  lag  network  for  which  the 

factor  0.0117  is  inherent  in  the  R-C  filter  design. 

Although  treatment  will  not  be  presented  here,  the 

algebraic  application  of  the  parameter  plane  technique  can 

be  readily  applied  to  combination  cascade  and  feedback 

compensat  ion . 

B.   DOMINANCY  OF  THE  SPECIFIED  ROOTS 

In  the   preceding  section  nothing  was  done  in  the 
calculations  to  make  the  specified  roots  a  dominant  pair. 
As  mentioned  earlier,  the  ability  to  predict  a  system's 
response  on  the  basis  of  the  location  of  a  pair  of  complex 
conjugate  roots  was  based  on  the  assumption  that  the 
magnitude  of  the  real  part  of  the  specified  roots  was  much 
less  than  that  of  the  remaining  roots  of  the  characteristic 
equation.   In  practice,  if  the  real  part  of  the  specified  or 
primary  roots  is  one  half  to  one  fifth  or  less  of  the  real 
parts  of  all  secondary  roots,  the  system  is  said  to  be 
dominant  in  the  primary  roots.   In  many  cases  the  system 
will  still  meet  the  specifications  even  if  two  pairs  of  complex 
roots  have  the  same  real  part,  provided  the  zetas  for  both 
pairs  of  roots  meet  the  specifications,  and  the  undamped 
natural  frequencies  are  such  that  the  component  time  responses 
are  not  highly  additive.   Further,  even  if  there  exists  a 
characteristic  root  whose  real  part  is  closer  to  the  origin 
than  that  of  the  primary  pair,  the  presence  of  a  closed-loop 
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zero  could  make  the  residue  of  the  close-in  root  negligible 
as  compared  to  the  primary  root.   If  possible,  however,  one 
usually  attempts  to  make  the  real  parts  of  all  secondary 
roots  large  in  magnitude. 

In  the   preceding  examples  it  should  be  pointed  out  that 
in  many  cases  there  were  actually  three  and  possibly  four 
variable  parameters.   For  instance,  the  forward  path  gain  was 
usually  a  fixed  value  in  the  computations  so  as  to  meet  the 
minimum  steady-state  accuracy  requirements.   There  is, 
however,  no  reason  why  the  gain  cannot  be  raised  above  the 
minimum  value,  thus  permitting  a  third  degree  of  freedom. 
When  cascade  and  feedback  compensation  are  used  simultaneously, 
the  forward  path  gain  and  tachometer  gain  become   the  third 
and  .fourth  parameters. 

Recall  that  the  system  characteristic  equation  has  the 

m     k 
form  f(S)  =   2   a.  S  =0  ,  where  a ,  =b,  a+c,  /3+d.  .   To  realize  the 

=0 

system  specifications,  one  places  a  pair  of  complex  roots 


at  S=-t  <u  +j<u -Vl-5?   ,  which  implies  that  S2+2  5  co  S+co2=0. 
Since  C   and  &)   are  known,  the  quadratic  can  be  divided 
out  of  the  characteristic  equation,  leaving  a  polynomial 
which  contains  the  secondary  roots  of  the  characteristic 
equation.   Since  only  two  of  the  variable  parameters  were 
used  in  fixing  the  primary  roots,  the  remaining  parameters 
will  appear  in  the  coefficients  of  the  quotient  polynomial, 
and  it  is  these  parameters  that  can  be  varied  to  achieve 
dominance. 
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Instead  of  division  to  obtain  the  quotient  polynomial, 
coefficients  of  like  powers  will  be  equated  to  achieve  a 
system  of  equations.   Let  the  quotient  polynomial  be 
given  by: 


n     k 
f,(s)  =  I       f,SK  =  0  (3-22) 

1      k=0   K 


where  n=m-2 ,  i.e.,  equation  (3-22)  is  of  order  two  less  than 
the  characteristic  equation.   Applying  equations  (2-1), 
(3-22),  and  the  quadratic  it  follows  that: 


2  2    n      k      m      k 

(S^+2U  S+u;)(  i       f,SK)  =  l      a,SK  (3-23) 

11    L      k=0   K       k=0   K 


Taking    a   =1    and    equating    coefficients   of    like   power: 


a      =    f    =1 

m  n 


am-l   "    fn-l+2ClUl 


2 
m-2  n-2      * l    l      n-1      1 


(3-24) 


;2  =  V2Wi+Vi 


2 
a      =2^^f+fjj 
1  1    1   o      1    1 


a0    =    ftfl 
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Equations  (3-23)  and  (3-24)  can  be  solved  for  the  coefficients 
f  in  terms  of  the  coefficients  a.   The  results  will  be 
applied  to  the  following  cases: 

Case  of  k=3,  n=l 
Equation  (3-23)  becomes: 


(S2+2C^1s+^)(f1S+f0)  =  S3+a2S2  +  a1S+a() 


Equating  coefficients  of  like  power  one  obtains: 


a3  =   f1  =  1 


a2    =    V2*lVl 


ai  =  V?+2CiVo 

f      2 

a0    -    fQ^ 

Solving  for  the  coefficients  f  results  in: 

fi  ■  1 

ao   (ai"'V 

f0  =  J  =  217^  ■  a2-2Cl"l  (3"25) 

1 


41 


Case  of  k=4,  n=2 
Proceeding  as  before,  the  a  coefficients  become: 


a4  -  f2  -  1 


a3  =  2^  +  fj 


a2  =  ^C^f^ 


al    =    flul+2Vlf0 


a   =  f  u? 
0    o  1 


(3-26) 


When  solved  for  the  coefficients  f,  equations  (3-26)  yield 


f2  =  X 


a2  Wl  a0  1  2^la0 

i1         a3    ^    x         2?1uj1         2£  240Ji3        co2      2      ^l 

+  co2    =    a0-2C  n^a0-^?+4C?u]2 


0         ^2         2C1^1         2CX  1  2       "1    1"3      1      "1    1 


(3-27) 
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Case   of    k=5 ,    n=3 
Similarly: 


a5    =    f3    =    X 


a      =    2£  oj    +f 
4  ^1    1      2 


2 
a      =0)       +2t  a     f+f 
a3  1  1    1X2    Xl 


a2=^l    V2*lWlfl  +  f0 


ai  =  "iV2  Wo 


ao  =  ^f0 


and: 


f3   "    l 


,2 


a9  1  T  1  ai 

f2    "    a4"2^l    "  ^2  +   a0(^   "  ^    "  ^3- 

1  111 

a3  al  +    a0         Wl 


2  ^1^1         2"Yl         ^1         2^ 
ai  ^  t  an  ?    9  9 

1 

f0        o)2 
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Although  the  coefficients  of  f  have  been  derived  for  only  up  to 
the  fifth  order  case,  they  can  easily  be  obtained  for  higher 
order  cases  if  necessary. 

Example  3-6  (Third  Order  Characteristic  Equation) 


±GhM 


T 


(s+p/  ) 

S+P 


S(S+4) 


Figure  3-4 

Design  a  cascade  compensator  for  the  system  of  Figure  (3-4) 
to  obtain: 

1.  Characteristic  roots  at  zeta=0.5  and  omega=40. 

2.  K  >250. 

e— 

3.  The  specified  roots  are  to  be  dominant. 
The  characteristic  equation  of  Figure  (3-4)  is: 


S3  +  (4+P)S2  +  (4P+Kr)S  +  KP  =  0 


or 


S3  +  (4+<*)S2  +  (4«+K/S)S  +K<*  =  0 


where  a=P  and  /?  =  y. 
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Here  G  =  ifi7  =  "277"'  so  eo=0>  ei=4>  e2=1'  V1'  V° 

Application  of  equations  (3-20)  yields: 


B   =  -K+co   =  -K+1600  B   =  -4cc  +  a)  2U   =  1440 


C  =    0  C2  =  ^jjK  =  -40K 


Dx  =  4co2  -  oj3U2  =  -57600       T>2    =  4oa  2U   -  a>  3U   =  6400 


From  equations ( 3-21)  are  obtained: 


a         57600 


-K+1600 


8         1440( -57600) -(-K+1600) (6400)  .  . 

-40K(-K+1600)  ^-^»; 


For  any  value  of  K,  equations  (3-28)  will  produce  values  of 
alpha  and  beta  that  provide  characteristic  roots  at  zeta=0.5 
and  omega=40.   However,  only  certain  ranges  of  K  will  meet 
the  steady-state  error  and  dominance  requirements .   To  satisfy 
the  error  consideration  it  is  necessary  that  K>1000.   Since 
f,(s)  is  of  order  one,  equations  (3-25)  apply  and: 


fl  -  X 

t0         ^2  '    2|  a)      a2  "  <4*la)l 
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From  the  characteristic  equation  it  is  seen  that 


a2  =  4+a 


a   =  4a  +  K 


a0  =  Ka 


The  real  part  of  the  specified  roots  is  |co=20.   Arbitrarily 
choosing  a  dominance  factor  of  five,  the  dominance  criterion 

becomes:   f~  >  5|cu  =  100.   To  satisfy  this  requirement,  the 

a0 
simplest  form  of  f~  will  be  chosen,  namely  f  =—77.   It  is 

0  o  ^2 

,,  .  Ka       57600K  1   36K     nn 

xnen  seen  tnat  iq  •  l600    i600( -K+1600)  1600-K 

This  requires  that  K>1176.5.   Since  K>1176.5  also  satisfies 

the  error  specification,  a  value  of  K=1180  is  arbitrarily 

chosen.   Using  this  value  of  K,  it  is  found  that: 

a=l37,  6=4.3  and  f  =101.   As  a  check,  from  the  expression 

f0=a2-2foo: 


fQ  =  (4+137)-2(0.5)(40)  =•  101 
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Example  3-7  (Fourth  Order  Characteristic  Equation) 


A 


(S+0.5)(S+1)(S  +  5)(S  +  10) 


T~ 1 


Figure  3-5 

Compensate  the  system  of  Figure  (3-5)  using  tachometer  plus 
acceleration  feedback  to  obtain: 

1.  Characteristic  roots  at  zeta  =  0.5  and  omega  =  2. 

2.  K  >12. 

e— 

3.  Dominance  of  specified  roots. 

The  characteristic  equation  of  Figure  (3-5)  is: 


S4+16.5S3+(73+a)S2+(82.5+B)S+25+K    =    0 


where    ^=KK      and    0=KK    . 
a  t 


G    = 


K 


S4+16.5S3+73S2+82.5S+25 
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By  inspection  eQ=25  ,  e  =82.5,  e2=73,  e3=16.5,  e.  =  l,  U  =1, 

U  =0,  and  U  =-1.   When  equations  (3-9)  are  employed  one  finds 


K-135        g  _  K-24 
4  2 


From  the  characteristic  equation  it  is  seen  that  a  =1,  a~=16.5 

a=73+<*,  a  =82.5+6,  a  =25+K.   Since  the  quotient  polynomial 

2 
f,(s)  is  a  quadratic,  i.e.,  S  +f  S+f=0,  equations  (3-27)  apply 

It  would  be  desirable  that ,  from  a  dominancy  viewpoint , 

f  -,  >5 1  ,^1  =  5.   However,  looking  at  the  dominancy  equations 

for  this  case  (equations  (3-27)),  it  is  seen  that  one  of 

the  several  expressions  for  f,  is  f  =a3-2|  w  ,  which  is  a 

fixed  constant  even  though  the  remaining  expressions  for  f 

involve   one  or  more  variables.   Thus,  f  =14.5.   Noting  the 

a0 
most  simple  expression  for  fn  in  equations  (3-27),  f  =  — ~. 

Now,  since  f  =14.5>5,  a  dominant  situation  already  exists. 

However,  the  system's  performance  can  be  further  improved 

by  choosing  appropriate  values  of  zeta  and  omega  for  the 

secondary  roots.   From  the  error  specification  it  is 

necessary  that  |-r>12  or  K>300.   Now  f  ( s )=S2+14 . 5S+-§ 

or  f  (s)=S2  +  14.5S+6.25  +  0.25K.   For  K=300,  tA&)  = 

S2  +  14.5S+81.25.   Therefore,  2I2oj2  =  14.5,  u  =81.25  implying 

oo  =9.   Then  zeta=0.806.   These  appear  to  be  reasonable 

values  for  £„  and  oo   since  the  secondary  roots  taken  aJone 

would  produce  less  overshoot  and  a  smaller  settling  time 

than  the  primary  roots.   Using  this  value  of  K  one  obtains 

for  alpha  and  beta: 
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a   =    41.25,  6   =    138 


and   since   a=KK.     and    6=KK    ,    IT     is    found   to    be    0.1375   while 
t  at 

K   =0.46. 
a 

As  an  added  bonus  of  this  method,  all  the  roots  of  the 
characteristic  equation  are  now  known  and  the  system's 
time  response  could  be  computed  if  desired. 
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IV.   PROGRAM  DESCRIPTION 

The  goal  in  developing  any  computer  aided  design  program 
should  be   twofold:   (1)  provide  the  user  with  a  single, 
easily  understandable,  easily  usable  and  comprehensive 
engineering  tool,  and  (2)  dramatize  its  efficiency  above 
that  of  other  currently  available  methods.   Through  the  examples 
of  the  following  chapter  the  second  goal  will  be  demonstrated. 
It  is  first  desirable  to  reveal  the  methodology  and  internal 
structure  of  the  parameter  plane  curve  program — as  a 
consequence  it  is  hoped  that  the  first  goal  will  be  affirmed. 

A.   THE  PROGRAM 

The  parameter  plane  curve--generating  program,  or 
"program"  as  it  will  be  called  henceforth,  consists  of  a 
large  driving  routine  which  includes  all  necessary 
calculations  with  which  to  generate  the  curve  data,  and 
several  supporting  subroutines  (i.e.,  curve  plotting,  root 
solving,  data  saving,  etc.).   This  entire  package  is 
included  as  a  user-selected  option  within  another  controls 
system  computer  aided  design  package.   Among  other  options, 
the  latter  CAD  program  includes  a  root-locus  analysis-- 
as  mentioned  earlier,  the  usefulness  of  either  the  parameter 
plane  or  root  locus  technique  for  design  of  a  controls 
system  is  somewhat  limited,  but  in  combination  their 
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effectiveness  is  synergistic  (Chapter  V  will  assert  the  dual 
roles  of  the  root-locus  and  parameter  plane  methods). 
Facilities  available  within  the  program  are  many;  the 
major  options  are: 

1.  Plotting  of  constant  zeta  curves,  with  zeta  as  a 
function  of  omega. 

2.  Plotting  of  constant  omega  curves,  with  omega  as  a 
function  of  zeta. 

3.  Plotting  of  constant  sigma  (real  root)  curves. 

4.  Plotting  of  constant  zeta-omega  curves. 

5.  Tabular  output. 

6.  Rescaling  of  the  plots. 

Input  of  certain  data  is  required  to  enable  the  program. 
These  inputs  include: 

1.  Starting  value  of  oo  . 

n 

2.  Decades  of  co   to  be  considered. 

n 

3.  Number  (and  values)  of  constant  zeta,  omega,  sigma, 
and  zeta-omega  curves. 

4.  Coefficients  associated  with  the  constant,  alpha,  and 
beta  terms  of  the  characteristic  equation. 

Each  of  the  basic  program  option  areas  will  be  described  in 

appropriate  detail,  as  well  as  their  interaction  with  the 

input  data. 

1.   Constant  Zeta  Contours 

In  practice  design  specifications  for  control 

systems  are  given  in  terms  of  percent  overshoot,  settling 

time,  error  constraints,  etc.   A  value  of  zeta  can  be 

associated  with  the  first  of  these  specif icat ions--that  is, 
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a  given  percent  overshoot  requirement  can  be  related  to  a 
specific  value  of  zeta.   Given  a  specific  zeta  value,  the 
program  calculates  the  alpha  and  beta  coefficients  of  the 
characteristic  equation  by  holding  the  zeta  value  constant 
and  varying  the  value  of  omega.   The  limits  within  which 
omega  is  varied  are  defined  by  the  user's  choice  of  the 
initial  <^  value,  and  the  number  of  decades  of  omega  to 
be  considered. 

From  the  nature  of  the  mapping  process,  it  is  clear 
that  when  the  contour  of  the  coefficient  plane  passes  through 
a  designated  point  (M-point),  the  original  mapping  contour 
on  the  S-plane  passes  through  a  point  which  is  a  root  of  the 
characteristic  equation.   The  zeta  value  chosen  for  the 
contour  is  then  the  zeta  for  the  root.   The  value  of  omega 
associated  with  the  M-point  is" the  radial  distance  from  the 
origin  of  the  S-plane  to  the  root.   Thus,  a  complex  root 
is  determined  when  the  M-point  lies  on  a  constant  zeta  curve 
of  the  parameter  plane.   The  value  of  this  root  and  its 
complex  conjugate  is: 


2 

S  =  -£C0+jO)/  l-$ 


If  the  characteristic  equation  is  such  that  several  complex 
roots  exist,  then  the  parameter  plane  curves  required  to 
realize  these  roots  must  all  pass  through  the  M-point.   If 
the  complex  roots  have  the  same  zeta  but  different  omega 
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values,  then  the  constant  zeta  curve  must  pass  through  the 
M-point  more  than  once.   In  fact,  once  any  point  on  the 
zeta  contour  is  defined,  omega,  alpha,  and  beta  are  also 
defined,  and  all  roots  of  the  characteristic  equation  are 
thus  fixed. 

2 .  Constant  Omega  Contours 

Within  the  program,  for  a  given  value  of  omega,  zeta 
is  varied  between  zero  and  one  inclusively  while  omega  is 
held  constant. 

As  with  the  constant  zeta  curves,  any  point  on  a 
constant  omega  contour  is  the  omega  for  a  complex  root  of  the 
characteristic  equation.   By  selecting  an  operating  point, 
zeta  is  also  defined  whereby  a  pair  of  complex  conjugate  roots 
is  established.   Again,  once  any  point  is  chosen  on  either 
a  constant  zeta  or  a  constant  omega  curve,  all  roots  of  the 
characteristic  equation  are  established. 

3 .  Constant  Sigma  Contours 

When  real  roots  are  to  be  evaluated,  it  is 
convenient  to  return  to  the  characteristic  equation: 

m     k 

X      a,SK  =  0  (2-1) 

k=0   K 

where,  again,  a,  represents  a  linear  combination  of  constant, 
alpha,  and  beta  terms.   S=-a  (a  real  number)  is  then  an 
equation  of  a  straight  line  on  the  parameter  plane.   If  any 
line  of  constant  sigma  value  passes  through  the  M-point, 


53 


then  the  alpha  and  beta  coordinates  of  the  M-point  satisfy 
the  characteristic  equation  for  a  real  root  located  at  -<r. 

For  program  considerations,  one  enters  a  positive 
value  for  sigma  (corresponding  to  a  real  root  at  -<0  and 
the  constant  sigma  contour  (straight  line)  is  developed. 
The  coordinates  of  any  point  on  this  curve  produce  a  real 
root  at  -a.      That  this  is  a  useful  tool,  consider  that 
system  specifications  can  be  achieved  by  placing  a  pair  of 
complex  conjugate  roots  of  the  characteristic  equation  at 
a  specific  location.   To  ensure  dominance  of  this  root 
pair,  the  real  part  of  the  complex  roots  so  placed  should 
be  smaller  in  magnitude  than  that  of  the  remaining  system 
roots.   Roots  placed  at  a  specified  sigma  value  can  thus 
ensure  at  least  one  real  root  whose  magnitude  is  greater 
than  the  real  part  of  the  intended  complex  conjugate  pair. 
4.   Constant  Zeta-Omega  Curves 

For  a  fixed  value  of  zeta  and  omega  a  pair  of  complex 
conjugate  roots  is  defined  in  terms  of  the  expression: 


2 

S  =  -£oj  +  joj/  l-£ 


The  real  part  of  these  roots  is,  thus,  defined  by  the 

zeta-omega  product.   Note  that  settling  time  is  defined  as 

4 
T   =  - —  .   If  the  £oo  product  is  known,  so,  too,  is  the 

S     COO 

duration  of  the  transient  response.   Thus,  by  specifying  a 
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constant  value  for  the  zeta-omega  product,  any  point  on  the 
contour  generated  by  this  value  will  produce  a  given 
settling  time. 

For  any  of  the  parameter  plane  contours,  it  is  desirable 
to  ascertain  the  values  of  the  characteristic  roots  for  every 
few  values  of  alpha  and  beta.   This  feature  has  been 
incorporated  within  the  tabular  output  facility  as  described 
below. 

5.   Tabular  Output 

For  each  of  the  zeta,  omega,  and  zeta-omega  parameter 
plane  curves,  an  arbitrary   though  reasonable  300  points  are 
calculated  with  which  to  plot  the  contours.   For  the  constant 
sigma  curves,  only  two  points  are  needed  to  define  the 
required  straight  lines  (in  practice,  4  points  are  generated 
to  ensure  that  the  sigma  contours  can  be  plotted  within  the 
user-defined  axes  limits).   Because  of  the  bulk  of  data 
points  so  generated,  tabular  output  is  offered  as  an  option 
(as  is  graphical  output),  and  all  points  so  generated  are 
listed  for  the  user.   In  addition,  it  is  worthwhile  to 
calculate  system  roots  for  given  values  of  alpha  and  beta. 
However,  computation  of  roots  for  each  alpha  and  beta  pair 
would  cost  unnecessary  computer  time  and  will  likely  tax  the 
user's  patience  with  the  bulk  of  output  so  generated.   Thus 
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the  characteristic  roots  are  generated  for  every  tenth 
pair  of  alpha  and  beta  values. 
6.   Plot  Rescaling 

Regardless  of  the  plot  scale  selected  by  the  user, 
the  program  generates  the  full  300  data  points  for  each 
curve  requested  (4  points  for  constant  sigma  lines).   When 
the  graphical  output  option  is  requested,  the  first  family 
of  curves  is  automatically  scaled  to  encompass  each  and 
every  data  point.   The  disadvantage  of  this  technique  is 
that,  because  most  activity  for  the  vast  majority  of  systems 
occurs  near  the  physical  origin  (i.e.,  alpha=beta=0) ,  the 
curves  may  at  first  appear  within  only  a  very  small  sector 
of  the  entire  plot  area,  and  often  they  are  indistinguishable 
from  one  another.   The  advantages  of  automatic  scaling  for 
the  first  set  of  parameter  plane  curves  far  outweigh  this 
disadvantage.   First,  by  plotting  all  available  data  points, 
the  possible  limits  for  alpha  and  beta  are  exposed--this  is 
important  if  very  large  values  of  alpha  and/or  beta  are 
required  to  meet  the  design  specifications.   Second,  for 
some  systems  the  area  of  activity  may  not  occur  near  the 
origin,  and  automatic  scaling  spares  the  user  the  task  of 
selecting  a  sector  and  possibly  missing  a  sector  of  interest. 


Although  a  seemingly  arbitrary  choice,  the  generation 
of  characteristic  roots  for  every  tenth  alpha,  beta  pair 
produces  a  very  tidy  output  on  the  common ly-used  IBM-3278 
computer  terminal. 
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Once  the  user  is  able  to  view  the  panoramic  alpha, 
beta  parameter  curves,  it  becomes  obvious  which  sector(s) 
are  of  interest.   The  user  then  has  the  option  to  rescale 
the  set  of  curves  by  selecting  upper  and  lower  limits  for 
the  alpha  and  beta  axes.   He  may  continue  to  rescale  the 
family  of  parameter  curves  as  often  as  is  desirable,  and  at 
any  time  the  autoscaling  option  may  be  recalled. 

The  curves  generated  by  the  above  program  are 
sufficient  to  explore  most  control  system  engineering 
problems.   The  use  of  these  parameter  plane  contours,  and 
their  interaction  with  one  another,  will  be  evidenced  in 
the  next  chapter.   The  source  code  listing  of  the  program 
is  included  as  Appendix  B. 

B.   INSTRUCTION  TO  THE  USER 

The  parameter  plane  program  is  highly  interactive — the 
user  is  prompted  for  each  required  input.   A  brief 
description  of  all  but  the  most  trivial  input  items  follows. 

-  Starting  value  of  w  :   For  most  control  systems  the 
initial  value  of  w  is  chosen  to  be  zero.   Because  & 
is  used  in  the  denominator  of  certain  of  the  parameter 
plane  equations,  ^  must  be  greater  than  zero.   However, 
the  user  may  choose  ^  arbitrarily  close  to  zero  if 
desired. 

-  Number  of  decades:   For  the  majority  of  control  system 
problems  a  suitable  number  of  decades  to  be  considered 
might  be  two  or  three.   For  higher  order  systems,  it 
would  be  advisable  to  start  with  a  slightly  larger 
number  of  decades,  especially  if  the  initial  ^  value 
is  small.   For  subsequent  families  of  curves,  the 
number  of  decades  can  easily  be  changed. 
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-  Constant,  alpha,  and  beta  coefficient  values: 
Characteristic  coefficients  are  requested  from  the 
highest  to  lowest  order  term.   By  way  of  an  example, 
a  third  order  characteristic  equation  might  be: 

S3  +  (3a  +  3  +  10)S2  +  as  +  (3+5)  =  0 

Here,  the  constant  coefficients  would  be  entered  in  the 
following  sequence:   1,10,0,5  while  alpha  and  beta 
coefficients  would  be  entered  as  0,3,1,0  and  0,1,0,1 
respectively . 

-  Zeta  values:   By  convention,    values  are  restricted  to 
between  zero  and  one,  inclusively. 

-  Sigma  values:   Positive  values  of  sigma  correspond  to 
negative  real  roots.   Since  few,  if  any,  practical 
engineering  applications  exist  for  designing  a  positive 
real  root  into  a  system,  negative  values  for  sigma  are 
disallowed. 

-  Omega  values:   Values  for  constant  omega  curves  are 
restricted  in  the  lower  limit  to  the  starting  w  value, 
and  in  the  upper  bound  by  go  xiodecades  .        n 

-  Zeta-omega  values:   As  with  the  constant  sigma  curves, 
values  for  constant  zeta-omega  contours  must  be  greater 
than  or  equal  to  zero. 

The  user  then  has  the  following  options: 

1.  Review  entries. 

2.  Change  any  entry. 

3.  Tabular  output. 

4.  Graphical  output. 

Remember  that  tabular  output  includes  300  data  points 
for  all  but  the  sigma  contours.   Characteristic  roots  are 
displayed  for  every  tenth  alpha,  beta  pair.   Because  of  the 
bulk  of  output  for  this  option,  use  it  only  when  necessary. 
If  a  printed  copy  of  the  tabular  output  is  desired,  type 
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"record  on"  before  invoking  the  program.   Upon  exiting  the 

program,  type  "record  off",  after  which  the  user  may  save 

the   preceding  terminal  session  in  a  listing  file  designated 

by  a  name  of  his  choosing.   Simply  print  the  listing  file, 

which  will  include  all  output  which  has  transpired  on  the 

terminal  between  the  two  calls  to  "record". 

When  graphical  output  is  requested,  all  curves  are 

superimposed  on  the  same  plot.   The  first  set  of  curves  is 

produced  with  an  autoscaling  feature,   which  plots  all 

points  calculated  (the  range  of  points  depends  on  your 

choice  of  initial  w  ,  number  of  decades,  zeta  values,  omega 

n  '  '     b 

values,  etc.).   For  most  characteristic  equations  only  the 
first  quadrant  (i.e.,  positive  alpha  and  beta  values)  will 
be  of  interest,  since  negative  values  usually  (but  not 
necessarily)  imply  negative  characteristic  coefficients  and, 
thus,  positive  roots  leading  to  system  instability.   The 
nature  of  the  first  (autoscaled)  set  of  curves  will  reveal 
the  actual  areas  of  interest  for  subsequent  plots  for  the 
same  system. 

Finally,  after  each  family  of  curves  is  plotted,  the 
user  has  six  additional  options: 

1.  New  problem. 

2.  Same  problem. 

3.  Root  finder. 

4.  Save  problem. 
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5.  Create  "DISSPLA  metafile". 

6.  Return  to  main  menu. 

Items  1,  4,  and  6  are  self-explanatory.   The  remaining 
options  deserve  some  additional  explanation. 

-  Option  2:   This  option  can  be  used  to  re-enter  the 
problem  at  a  point  prior  to  actual  plotting.   Then, 
specific   input  values  can  be  added  or  modified, 
tabular  and/or  graphical  output  can  be  requested,  and 
entries  can  be  reviewed.   It  is  within  this  option  that 
the  graph  coordinate  axes  can  be  rescaled.   If  user- 
defined  scaling  is  desired,  the  minimum  and  maximum 
values  for  the  axes  are  requested. 

-  Option  3:   Although  within  the  tabular  output  feature  a 
set  of  characteristic  roots  is  produced  for  every  tenth 
pair  of  alpha  and  beta  values,  the  bulk  of  output  using 
that  option  may  prove  excessive  for  some  applications. 
Here,  the  user  has  the  option  of  choosing  specific 
values  of  alpha  and  beta  (e.g.  extracted  from  the 
family  of  parameter  plane  curves)  and  obtaining  the 
system  roots. 

-  Option  5:   The  program  provides  a  choice  frbm  among 
four  graphic  output  devices.   Usually  the  user  will 
nominate  the  TEK618  graphics  terminal  due  to  its 
relatively  high  quality  plot  resolution.   Once  the 
user  has  produced  a  parameter  plane  plot  to  his  liking, 
he  may  wish  a  final  plot  of  very  high  resolution.   By 
selecting  this  option,  the  plot  is  stored  as  a  DISSPLA 
metafile,  and  the  program  is  terminated  (termination 

of  the  program  is  necessary  at  this  point  due  to  an 
anomaly  of  the  DISSPLA  graphics  package).  Simply  type 
"DISSPOP"  and  follow  the  simple  instructions,  choosing 
the  default  options  as  they  are  presented.  Within  the 
"DISSPOP"  routine,  any  of  several  output  devices  can  be 
called,  including  the  high  resolution  Versatec  plotter 
and  the  3800  laser  printer. 
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V.   PARAMETER  PLANE  CURVES-GRAPHICAL  METHOD 

A.   GRAPHICAL  SOLUTION 

The  algebraic  solutions  discussed  in  Chapter  III  have  the 
disadvantage  that  a  fixed  value  of  zeta  and  omega  must  first 
be  chosen  to  compute  the  alpha  and  beta  terms.   In  some 
instances  it  is  possible  to  modify  the  remainder  polynomial 
so  as  to  ensure  that  the  specified  roots  are  dominant. 
However,  it  is  not  always  possible  to  guarantee  that  roots 
placed  at  a  specified  location  can  be  made  dominant.   Thus, 
an  exhaustive  trial-and-error  procedure  may  be  required  to 
achieve  the  best  values  for  the  various  parameters.   Trial- 
and-error  may  also  be.  required  in  the  design  of  cascade 
compensators  where  a  specific  root  location  may  require 
parameter  values  that  are  not  physically  realizable.   In 
these  cases,  the  calculation  must  be  repeated  in  terms  of 
slightly  modified  specifications;  possibly  a  different  means 
of  compensation  must  be  used. 

To  avoid  this  trial-and-error  analytical  approach  ,  one 
can  employ  the  graphical  solution.   Once  a  family  of  curves 
is  generated  by  the  program  one  can,  by  choosing  an  M-point 
in  the  parameter  plane,  obtain  from  the  curves  the  n  roots 
of  the  nth  order  characteristic  equation.   The  trial-and-error 
procedure  can  then  be  done  visually  to  reveal  an  operating 
point  which  best  meets  the  given  specifications. 
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Example  5.1  (An  Attitude  Control  System  for  Large  Launch 
Vehicles) 
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Figure  5-1 

Figure  5-1  shows  the  mechanization  for  a  control  system  for  a 
large  launch  vehicle;  Figure  5-2  shows  the  equivalent  block 
diagram 


Figure  5-2 
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where  the  equivalent  G(s)  is: 


(aQS2  +  K  S  +  K2)(S2-C) 
G(s)  =  


3        2 
S3(S_  +  2£^S   +  g  + 

CO        CO  l 

e      e 


From  G(s)    one   obtains   the   system's    characteristic    equation 


q6  2?                                               . 

^2  +  oT^  S      +    (l+aQ)S*   +    (a^+K    )S°   +    (K2-aQC)SZ 

co  e 
e 


-   CK   S    -    CK2    =    0 


where   a~  ,  a, ,  K. ,  K„  =  control  system  gains 
K       =  damping  ratio   of  control  servo 


to   =  natural  frequency  of  control  servo 


Gains  K   and  K„  are  chosen  as  the  system  parameters  a  and  B  , 

respectively,  to  be  portrayed  in  the  parameter  plane.   One 

must  then  find  suitable  values  for  these  parameters  to  yield 

a  desired  stability  margin  and  a  satisfactory  transient 

response.   For  a  typical  choice  of  system  parameters  (such 

as  those  describing  Saturn  V),  it  is  assumed  that: 

K       =   0.717 
e 

w   =  4.71  Hz  =  29.594  radians 
e 

a„  =  -0.5 
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a   =  1.0 
C  =  0.7 


Then  the  system's  characteristic  equation  becomes: 

0.0011S6  +  0.0485S5  +  0.5S4  +  (0.7+ct)S3  +  (0.35  +  6  )S2 


-  0.7CS  -  0.73  =  0 

Various  values  of  a~  and  a1  were  used  to  deduce  their  effect 
on  the  stability  regions,  which  is  indicated  in  Figure  5-3. 
The  numbers  of  stable  and  unstable  roots,  respectively,  are 
portrayed  in  parentheses  for  each  region  of  the  parameter 
plane. 

The  analysis  procedes  as  follows:   the  £ =0  contour  repre- 
sents boundaries  of  stability  (or  relative  stability  when 
K> 0)    associated  with  pairs  of  complex  conjugate  roots.   The 
S=0  (sigma=0)  curve  represents  real  root  stability  boundaries 
The  region  of  stability  is  thus  that  area  bounded  by  these 
two  contours.   See  Figures  5-4a  and  5-4b  for  a  magnified 
view  of  this  area.   Any  negative  alpha-beta  pair  from  within 
this  region  will  exhibit  six  stable  (i.e.,  all  within  left- 
half  S-plane)  roots  and  system  stability  will  be  assured. 
Note  that  from  the  form  of  the  characteristic  equation, 
both  alpha  and  beta  must  be  negative  to  obtain  a  stable 
system.   To  illustrate,  let  us  select  an  arbitrary  operating 
point  constrained  to  lie  within  the  lower  loop  of  Figure  5-4b 
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Figure    5-3 
Parameter   Plane   Curves    for 

0.0011    S6    +    0.0485S5    +    0.5S4    +    (0.7+a)S3    +    (0.35+S)S2 

-    0.73 S    -    0.7a    =    0 
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ATTITUDE  CONTROL  SYSTEM 
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Figure    5-4a 
r.  Para.meter   Plane   Curves    for 

0.0011S      +    0.0485S5   +    0.5S4    +    (0.7+a)S3    +    (0.35+3)S^ 
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ATTITUDE  CONTROL  SYSTEM 
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Figure   5-4b 
Parameter   Plane   Curves    for  o  o 
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and  also  satisfying  the  requirement  that  both  alpha  and  beta 
be  negative.   Superimposed  on  Figure  5-4b  is  the  constant 
C=0.5  contour,  upon  which  our  M-point  might  be  chosen.   If 
we  select,  say,  a=-0.15  and  B=-0 .02  ( corresponding  to  5=0.5 
and  w=0.5),  the  system  roots  are: 


-dominant  roots 


-0.26  +  jO.45 
-0.26  -  jO.45 
-0.32  +  jO.12 
-0.32  -  jO.12 
-14.6  +  jO 
-26.7  +  jO 

Coincidentally ,  the  roots  associated  with  our  choice  of  zeta 
and  omega  are  seen  to  be  dominant.   The  actual  choice  of  an 
operating  point  may  depend  on  other  criteria  not  discussed 
here. 

Since  K   and  K   (a  and  3  ,  respectively)  are  functions 
of  w  ,  £  ,  a    a   C,  &  f    and  E,     ,  they  may  be  determined  for 
various  instances  of  flight  by  plotting  several  constant 
zeta  and  constant  omega  curves.   Actually,  K,  and  K„  vary  so 
little  within  the  range  of  values  used  for  C,  for  specified 
values  of  5  and  '■&>    that  it  becomes  possible  to  choose  constant 
values  for  K   and  K  . 

This  has  been  a  relatively  simplistic  treatment  of  a 
complicated  control  system  problem,  but  it  demonstrates  the 
power  of  the  parameter  plane  graphical  technique.   Knowing 
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nothing  more  than  the  system's  characteristic  equation,  the 
static  and  dynamic  stability  boundaries  may  be  obtained  to 
define  the  area  of  overall  stability.   Of  course,  other 
system  constraints  may  exist  to  further  limit  this  area. 
As  a  note,  the  root-finding  option  available  within  the 
program  was  used  to  confirm  the  numbers  of  stable  and 
unstable  roots  lying  within  each  region  of  Figures  5-3  and  5-4 

Example  5-2  (Alternator  Voltage  Regulator) 
In  designing  a  voltage  regulator  of  the  type  shown  in  Figure 
5-5,  we  must  find  values  for  K  ,  K  ,  and  K9  that  provide 
stability  and  good  transient  performance. 


^ 


EXCITER 


K2S    -HK^+Kq 
(T.S-HXT-S-H) 


ALTERNATOR 


Figure  5-5 

The  characteristic  equation  of  the  system,  including  alternator, 
tie-line,  etc.  is  (see  reference  8  for  details): 

0.0095S5  +  0.1325S4  +  1.72S3  +  (K  +7.55)S2  +  (S  +9,1)S 

«  J. 

+  KQ-2.5  =  0 
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We  will  consider  K  =0  as  fixed  and  a   and  3  as  variable  para- 

o  * 

meters  corresponding  to  K  +7.55  and  K  +9.1,  respectively. 

Curves  of  a  versus  3  are  plotted  in  Figure  5-6  with  oj   as 

*  to  n 

the  variable  parameter  (for  this  system,  K   and  K_  are 
known  to  be  functions  of  co   alone).   Since  the  program  plots 
constant  zeta  curves  as  a  function  of  varying  omega,  these 
curves  were  selected  as  logical  candidates  to  study  the 
problem. 

The  beta  axis  corresponds  to  the  zero  value  of  damping 
constant  (sigma)  since  below  the  zeta=1.0  curve,  the  real 
roots  (Thaler  and  Brown  1960)  are  the  negative  slopes  of  the 
tangents  drawn  from  the  point  in  question  to  the  £=1  curve. 
The  machine  is  then  stable  for  any  (a,  3)  pair  between  the 
C=0  and  S=0  curves.   The  greatest  stability  of  the  machine 
is  then  possible  when  both  zeta  and  sigma  are  largest. 
Further,  the  best  stability  can  be  expected  in  the  region 
bounded  by  the  £=0.3  and  £=1.0  contours  (Kabriel  1967). 
Similar  stability  limits  of  a  and  3  can  be  investigated  by 
taking  K   as  10,  20,  30  and  so  on. 

Consider  now  K  =0  as  fixed.   Then  a  and  3  will  represent 
the  variable  parameters  K  +7.55  and  K^-2.5,  respectively. 
The  characteristic  equation  then  becomes: 

0.0095S5  +  0.1325S4  +  1.72S3  +  aS2  +  9. IS  +3=0 


70 


A 


gain,  K„  and  K„  are  known  to  be  functions  of  ui   alone.   Here 
°       0       2  n 


however,  co   can  be  solved  explicitly  [Ref.  8]  and  one  obtains 

three  straight-line  equations: 

3=0 
5.421a    -3    =    3.894 
175.63a    -3    =    4085.0 

or  in  terms  of  system  parameters: 


KQ  =  2.5 

5.421(K0+7.55)-K  =  3.894  -  2.5 
2         o 

175.63(K  +7.55)-KQ  =  4085.0   -  2.5 


These  are  plotted  in  Figure  5-7.   The  triangular  region 
bounded  by  these  three  lines  represents  a  stable  region.   The 
values  of  (a ,  3 )  within  the  triangle  and  hence  corresponding 
values  of  K~  and  K„  can  be  predicted  for  stable  operation. 
The  procedure  can  be  repeated  for  further  investigation  by 
taking  K  =15,  30,  45,  etc. 

The  analytical  parameter  plane  technique  can  be  used  to 
determine  the  stability  limits  of  K„  and  K   by  choosing  several 
constant  values  of  K  .   But  since  K   has  to  be  selected 
arbitrarily  for  this  purpose,  this  method  becomes  cumbersome 
and  time-consuming.   The  method  presented  here,  on  the  other 
hand,  can  be  used  to  determine  the  stable  range  of  either 
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ALTERNATOR  VOLTAGE  REGULATOR 
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Figure  5-7 
c-     Parameter  Plane  Curves  for 
0.0095S   +  0.1325S4  +  1.72S3  +  aS2  +  9. IS  +3=0 
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K   and  K9  or  K   and  K   fixing  the  third  parameter.   Over  and 
above  the  ability  to  predict  stable  operation,  the  method 
povides  a  direct  measure  of  the  damping  at  and  around  a 
chosen  operating  point. 

Example  5-3  (Lead  Compensator) 
Consider  the  system: 


J®, 


K(S-t-Z) 


S+P 


1 


S(S-H)(S-»-4) 


The  characteristic  equation  is: 


S4  +  (5+P)S3  +  (4+5P)S2  +  (4P+K)S  +  KZ  =  0 


We  choose  to  cancel  the  pole  at  S=-l  with  the  zero;  thus 
Z=1.0  and  the  characteristic  equation  becomes: 

S4  +  (5+P)S3  +  (4+5P)S2  +  (4P+K)S  +  K  =  0 

Let  P  =  oc  and  K  =  3.   Then  the  parameter  plane  curves  are  as 
shown  in  Figure  5.8,   For  the  coefficients  of  the  character- 
istic equation  to  remain  positive  (and  thus  ensure  stability), 
it  is  convenient  to  consider  only  positive  values  for  alpha 

and  beta.   By  inspection  we  choose  ^  =  0.5  and  cj  =2.0  as  a 

'     *  n 
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Figure  5-8 
Parameter  Plane  Curves  for 

S4  +  (5+ct)S3  +  (4+5ct)S2  +  (4a+3)S  +  3=0 
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LEAD  COMPENSATOR 


Ld 
CD 


2.5      5 

ALPHA 


Figure  5-9 
Parameter  Plane  Curves  for 

S4  +  (5+a)S3  +  (4+5a)S2  +  (4^+3  )s+ 3 


=  0 
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good  operating  point,  for  which  a  =  4.0  and  3  =24.0.   The 
corresponding  vicinity  of  Figure  5-8  is  re-scaled  in  Figure 
5-9.   Again,  using  the  root-finding  option,  the  roots 

associated  with  this  (?  ,  w  )  pair  are  shown  to  be  dominant. 

'   n 

Example  5-4  (Lag  Compensator) 
If  we  are  especially  concerned  with  steady-state  accuracy 
for  a  ramp  input,  it  may  be  advisable  to  design  a  lag 
compensator.   The  parameter  plane  permits  us  to  consider 
steady-state  accuracy  while  designing  the  transient  response. 
If  our  system  is  the  same  as  that  considered  for  the  lead 
compensator,  the  characteristic  equation  remains: 

S4  +  (5+P)S3  +  (4+5P)S2  +  (4P+K)S  +KZ  =  0 

KZ 
But  now  the  error  coefficient  is  K   =  -r=r    .   Having  three 

e    4P 

unknown  parameters,  K,  Z,  and  P,  two  must  be  selected  (or 
some  combination  of  two)  to  be  ct  and  3  while  a  numerical  value 
is  assigned  to  the  third.   Once  this  choice  has  been  made 
the  parameter  plane  curves  can  be  calculated  and  plotted,  and 
the  loci  of  constant  K   can  be  superimposed.   Let  Z=0.1, 
P= a,  and  K= 3.   The  characteristic  equation  becomes: 

S4  +  (5+ot)S3  +  (4+5ct)S2  +  (4a+3)S  +  0.13  =  0 


Parameter  plane  curves  are  shown  in  Figures  5-10  and  5-11. 

0  1& 

Lines  of  K   =   ' a  =  0.1,  0.2,  0.3,  etc.  may  be  superimposed. 

If  we  select  K   =0.15  and  (similarly  to  the  lead  compensator 
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Figure    5-10 
.  Parameter   Plane   Curves    for 

S      +    (5+a)S3    +    (4+5a)S2    +    (4a+6)S    +    O.lB    =    0 
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LAG  COMPENSATOR 
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Figure    5-11 
4  Parameter   Plane   Curves    for 

S      +    (5+a)S3    +    (4  +  5a)S2    +    (4a+B)S    +    0.13    =    0 
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example)  C=0.5  and  00  =2.0,  then  a=2.14  while  3  =12  .  89 .   This 

produces  dominant  roots  at  -1.0  +  jl.7. 

Let  us  reconsider  the  systems  for  which  a  lead  compensator 

and  a  lag  compensator  were  designed.   The  characteristic 

equation  as  well  as  the  error  coefficient  each  contain  three 

unknowns  (if  we  assume  a  numerical  value  for  K  ).   We  can 

e 

imbed  the  error  coefficient  in  the  characteristic  equation 

by  direct  substitution,  thereby  eliminating  one  of  the 

unknowns.   Let  us  eliminate  the  gain  parameter  K  -  note  that 
4PK 

Returning  to  the  characteristic  equation: 

4  o  o        4PK 

S   +  (5+P)S   +  (4+5P)S   +  (4P+-=-^)S  +  4PK   =  0 

&  e 


p 
Letting  P  =  a  and  =  =  3  ,  the  characteristic  equation  becomes 

S4  +  (5-wc  )S3  +  (4+5a  )S2  +  (4a+4K  3  )S  +  4K  a  =0 

e        e 


If  a  constant  value  is  now  chosen  for  K  ,  say  K   =  2.0,  the 

e  *         J      e       ' 

characteristic  equation  becomes: 

S4  +  (5-H*)S3  +  (4+5ct)S2  +  (4ct+83)S  +  8ot   =  0 


for  which  the  parameter  plane  contours  are  first  shown  in 
Figure  5-12  and  are  further  magnified  in  Figure  5-13.   Now 
when  any  operating  point  is  chosen  on  the  parameter  plane 
curve(s),  the  selected  (a,3)  pair  generates  K  =2.0.   Of 
course,  this  procedure  can  be  repeated  for  any  choice  of  K 
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EMBED  ERROR  COEFFICIENT   (KE=2) 
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Figure    5-12 
4  Parameter   Plane   Curves    for 

S      +    (5+a)S3    +    (4+5«)s2    +    (4cw-83)S    +    8a  =    0 
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EMBED  ERROR  COEFFICIENT   (KE=2) 
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Figure    5-13 
4  Parameter   Plane    Curves    for 

S      +    (5+a)S3    +    (4+5a)S2    +    (4a+8B)S    +    8a   =    o 


10 
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On  the  parameter  plane  curves  let  us  choose  =  =0.5  and 


U    =  2.0 
n. 

for   which: 

a    = 

6.00    =    P 

S    = 

4.50    =   | 

a 

1.33    =    Z 

K    = 

83         =    36 

KZ 


As  a  check,  Kg  =  -^   =  2.0.   The  roots  associated  with  the 
selected  zeta  and  omega  values  are  shown  to  be  dominant 
using  the  root-finding  algorithm  and  are: 


■1.00    +    jl.73 
■1.00    -    jl.72 
■1.63    +    jO 
•7.37    +    jO 


-  Dominant  Pair 


Example  5-5 
As  a  final  engineering  example,  consider  the  system  belo 


w 


K 


S(S+10)(S2+5S-t-100) 


KS 


> 
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For  this  system,  the  following  specifications  are  to  be  met: 

1.  Set  K  to  the  stability  limit. 

2.  Place  a  dominant  root  pair  within  the  following  region: 
0.4<?<_0.7,  and  2<_oj^6. 

3.  Both  tachometer  and  acceleration  feedback  may  be  used, 
but  if  possible  choose  only  one. 

From  the  figure  the  uncompensated  system's   open  loop  transfer 

function  is: 

GH  -  -1  -  — 


S(S+10)(S2+5S+100) 


From  which,  when  expanded,  the  characteristic  equation 
becomes: 

S4  +  15S3  +  150S2  +  1000S  +  K  =  0 

To  determine  the  value  of  K  at  the  stability  limit  the  Routh 
array  is  employed: 


1 

150 

K 

15 

1000 

0 

1250 

15K 

0 

1.25xl06-2 

25K 

0 

0 

15K 

0 

0 
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Here,  the  stability  limit  is  seen  to  be  K=5555.5.   If  both 
tachometer  and  acceleration  feedback  are  used  the  compensated 
system's  characteristic  equation  becomes: 

S4  +  15S3  +  (150+5555. 5«)S2  +  ( 1000+5555 . 53 )S  +  5555.5=0 

where  a=K  ,  $=K  ,  and  K  has  been  set  to  the  stability  limit. 
a     l 

Parameter  plane  curves  for  this  system  are  shown  in  Figure 

5-14. 

From  these  curves,  the  following  analysis  can  be  made. 

The  origin  of  the  parameter  plane  corresponds  to  the  roots 

of  the  uncompensated  system.   Since  the  5=0  curve  passes 

through  the  origin,  two  roots  are  located  on  the  ju>  axis  of 

the  S-plane  as  was  to  be  expected  from  the  Routh  array.   The 

remaining  two  roots  are  also  complex  and  correspond  to  ?=0.8 

and  co=5.0.   It  is  important  to  note  that  when  an  operating 

point  involves  two  different  pairs  of  complex  roots,  then  the 

curves  for  two  different  values  of  omega  and  two  different 

values  of  zeta  must  pass  through  the  point.   To  determine 

which  value  of  omega  corresponds  to  which  value  of  zeta,  it 

becomes  necessary  to  refer  to  the  program's  tabular  data 

output,  which  is  not  included  here  due  to  lack  of  space. 

With  K  =0,  the  effect  of  tachometer  feedback  alone 
a 

corresponds  to  movement  of  the  M  point  along  the  3  axis. 
In  Figure  5-14  the  unstable  region  is  determined  by  an 
inspection  of  the  way  the  constant  zeta  curves  tend 
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Figure  5-14 
Parameter  Plane  Curves  for 

S4  +  15S3  +  ( 150+5555. 5ct  )S2  +  ( 1000+5555 . 5  S)S  +  5555.5  =  0 
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as  zeta  increases.   Since  the  3  axis  is  always  in  the 
unstable  region,  it  is  concluded  that  tachometer  feedback 
alone  cannot  stabilize  the  system. 

The  effect  of  acceleration  feedback  alone  can  be  observed 

by  traveling  along  the  a-axis  of  Figure  5-14.   If  K   is 
varied  between  0.01  and  0.06,  the  system  will  exhibit  two 
pairs  of  complex  roots  with  the  following  ranges  of  values 
for  zeta: 

0.3<?<0.5    and    0.25<5<0.32 

If  both  tachometer  and  acceleration  feedback  are  used, 
it  is  seen  that  tachometer  feedback  will  in  general  cause 
the  zeta  of  one  pair  of  roots  to  increase  while  the  zeta 
of  the  other  pair  decreases.   Acceleration  feedback  alone 
would  appear  to  be  the  better  choice. 

From  the  set  of  curves  it  is  determined  that  with  K  =0 
and  K  =0.012,  four  complex  roots  are  located  with  associated 

ct 

values  at  £=0.45,  o>=4.0,  and  £=0.32,  oo=13.0.   Since 
(0.45)(4)=18<< (0.32)(13)=4.15,  it  is  apparent  that  the  roots 
at  5=0.45  and  ^=4.0  are  dominant.   The  specifications  have 
been  met  and  the  problem  is  solved. 

B.   MISCELLANEOUS  ASPECTS  OF  THE  PARAMETER  PLANE 

It  can  be  demonstrated  that  constant  zeta  parameter 

plane  curves  of  order  two  through  five  originate  at  a  point 

M  N 

where  alpha  =  r=  and  beta  =  ^-,  where  M,  N,  and  K  are 
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determined  by  the  zero  and  first  power  coefficients  only. 
Intuition  can  be  used  to  conclude  that  constant  zeta  curves 
of  any  order  originate  at  this  common  point  which  is 
determined  only  by  the  zero  and  first  power  coefficients. 
An  exception  is  when  K=0.   In  this  case  the  origin  point 
depends  on  higher  order  coefficients  and  its  location  will 
be  obvious  given  a  specific  problem.   If  K  is  not  zero, 
the  origin  is  independent  of  the  order  of  the  characteristic 
equation . 

Inspection  of  the  expressions  for  alpha  and  beta  indicates 
that  the  shape  of  the  constant  zeta  curves  as  omega  becomes 
larger  is  primarily  determined  by  the  coefficients  of  higher 
power,  and  in  general  the  curves  become  more  complex  and 
less  well  behaved  as  the  order  of  the  characteristic 
equation  increases.   For  a  given  characteristic  equation, 
an  increase  in  complexity  can  be  observed  as  alpha  and  beta 
appear  in  more  coefficients. 

All  constant  zeta  curves  tend  to  plus  or  minus  infinity. 
The  relative  magnitudes  of  the  coefficients  determine  whether 
the  limit  is  plus  or  minus  infinity.   It  is  therefore 
necessary  to  choose  a  frequency  range  of  interest  before 
plotting  the  curves,  thus  limiting  the  analysis  to  one 
"window"  of  the  infinite  plane. 

Since  no  stability  criteria,  either  relative  or 
absolute,  has  been  established  for  the  parameter  plane,  it 
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is  necessary  to  base  the  stability  analysis  on  observing 
which  way  the  curves  tend  as  omega  and  zeta  are  varied. 
For  this  reason  it  is  worthwhile  to  plot  curves  for  as  many 
values  of  zeta,  omega,  sigma,  and  if  desired,  zeta-omega, 
as  are  necessary  to  ascertain  the  pattern. 
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VI.   CONCLUSIONS 

Parameter  plane  techniques  have  been  applied  to  the 
compensation  of  linear  control  systems.   General  equations 
have  been  derived  for  the  cases  of  feedback,  cascade,  and 
combination  feedback-cascade  compensation,  to  enable  one  to 
place  a  pair  of  complex  conjugate  roots  at  a  specific 
location  in  the  S-plane,  while  simultaneously  satisfying 
the  steady-state  accuracy  requirements.   A  dominancy 
technique  has  been  introduced  whereby  once  a  pair  of  complex 
roots  is  fixed,  the  remaining  roots  of  the  characteristic 
equation  can  be  manipulated  to  ensure  that  the  specified 
roots  are  dominant. 

The  impetus  for  development  of  a  parameter  plane  program 
was  to  provide  the  user  with  a  quick,  simple  means  of 
obtaining  the  information  available  in  the  analytical 
solution  to  control  system  compensation,  while  avoiding  the 
painstaking  labor  of  trial-and-error  analysis  inherent  in 
that  technique.   Several,  practical  engineering  examples 
have  been  presented  to  demonstrate  the  superiority  of  the 
graphical  technique.   To  date,  no  other  package  is  known  to 
offer  the  fully  interactive  and  comprehensive  capabilities 
of  the  parameter  plane  program. 
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By  itself  the  program  allows  one  to  design  a  control 
system  compensation  model  for  most  systems.   However,  for 
some  lightly  damped  systems  containing  mechanical  resonances, 
the  amount  by  which  zeta  or  omega  are  incremented  in  the 
parameter  plane  equations  may  be  so  large  as  to  not  detect 
the  resonance  peaks.   This  information  would  be  available  from 
either  a  root-locus  or  Bode  analysis.   For  still  other 
systems,  one  might  be  interested  in  the  way  the  roots  of  the 
characteristic  equation  extend  from  the  open  loop  poles  and/or 
zeros.   Since  the  parameter  plane  equations  are  calculated 
using  only  the  characteristic  equation,  no  knowledge  of  open 
loop  poles  or  zeros  is  available.   Again,  a  root-locus  method 
would  reveal  this  information.   Incorporated  into  one 
comprehensive  package  which  includes  Bode  and  root-locus 
analyses,  the  program  provides  the  capability  to  investigate 
the  entire  gamut  of  linear  control  system  architecture. 

A  basis  for  further  investigation  involves  plotting  the 
parameter  plane  contours  for  systems  that  are  non-linear 
in  the  alpha  and  beta  terms--i.e.,  those  systems  which  contain 
alpha-beta  product  terms.   Although  the  recursion  technique 
used  in  this  text  has  distinct  advantages  over  the  matrix 
approach  for  the  linear  case,  as  pointed  out  earlier,  the 
matrix  technique  would  be  the  method  of  choice  for  the  non- 
linear case. 
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APPENDIX  B 
PARAMETER  PLANE  PROGRAM 


100 


101 


102 


103 

104 
105 

106 


SUBROUTINE  LPARAM 

DIMENSION  A(350),  B(350), 
X      AG(350),  BGC350),  AJ(IOO),  BJ(IOO),  CJ(IOO), 
X      ZETAUOO),  SIGMAC100),  W(lOO),  ZW(IOO) 

INITIAL  ASSIGNMENTS 
CHARACTER**   SCHAR/'S  »/,YES/'Y  '/,NOO/'N  '/,  BLANK/'   •/, 
X      ENCHAR/'E  '/,  WNCHAR/  •  WN'/  ,  NDCHAR/  'ND'/,  NOCHAR/  '  NO  '/, 
X      AJCHAR/'AJ'/,BJCHAR/'BJ'/,CJCHAR/'CJ'/,NSCHAR/'NS'/, 
X      NZCHAR/ ' NZ ' / , ZWCHAR/ ' ZW  '/ , NWCHAR/ • NW » / , PRCHAR/ ' PR • / , 
X      NCCHAR/'NC'/,TICHAR/'TI'/ 
CHARACTER**   TABLE,  GRAPH,  CHANGE,  REPLY,  OPT,  LABEL(9) 
COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  C J ,  AJ,  BJ, 
X  NZ,  ZETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

DATA  ENTRY  FROM  FILE  OR  CONSOLE? 
MINMAX  =  1 
GRD  =  0. 
CHANGE  =  BLANK 
CALL  EXCMSCCLRSCRN') 
WRITE(6,500) 
CALL  READC  (REPLY) 
IF  (REPLY  .NE.  "D*)  GOTO  101 

CALL  GETIT 

MINMAX  =  0 

GO  TO  200 
CONTINUE 

GET  CURVE  TITLE 
CONTINUE 

CALL  EXCMSCCLRSCRN') 
WRITE(6,501) 
CALL  READL  (LABEL) 
CALL  ASTER  ( LABEL , LABEL) 
IF  (  CHANGE  .EQ.  TICHAR  )  GO  TO  200 

GET  STARTING  VALUE  OF  WN 
CONTINUE 
WRITE(6,502) 
CALL  READR  (WN) 

IF  (WN)  104,104,105 

WRITE(6,503) 

GO  TO  103 
IF  (  CHANGE  .EQ.  WNCHAR)  GOTO  200 


107 


CONTINUE 
WRITE(6,504) 
CALL  READI  (ND) 
IF  (  CHANGE  .EQ 


CONTINUE 

WRITE(6,505) 

CALL  READI  (N02) 

NC  =  N02+1 

IF  (CHANGE  .EQ.  NOCHAR  )  GOTO  200 


GET  THE  NUMBER  OF  DECADES  CONSIDERED 

NDCHAR  )  GOTO  200 

GET  THE  ORDER  OF  THE  CHARACTERISTIC  EQN 
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c 

C  GET  THE  NUMBER  OF  CONSTANT  ZETA  CURVES 

108  CONTINUE 

CALL  EXCMS  CCLRSCRN') 

WRITE(6,506) 

CALL  READI  (NZ) 

IF  (NZ  .LT.  1)  GOTO  110 
C 
C  GET  THE  VALUES  OF  ZETA 

WRITE(6,507) 

DO  110  I  =  1,NZ 

109  WRITE(6, 508)1 

CALL  READR  (ZETA(I)) 

IF  (ZETA(I)  .LT.  "0.  .OR.  ZETA(I)  .GT.  1.)  WRITE(6,509) 

IF  (ZETA(I)  .LT.  0.  .OR.  ZETA(I)  .GT.  1.)  GO  TO  109 

110  CONTINUE 

IF  (  CHANGE. EQ.  NZCHAR  )  GOTO  200 
C 
C  GET  THE  NUMBER  OF  CONSTANT  SIGMA  CURVES 

111  CONTINUE 

CALL  EXCMS  CCLRSCRN') 

WRITE(6,510) 

CALL  READI  (NS) 

IF  (NS  .LT.  1)  GOTO  113 
C 
C  GET  THE  VALUES  OF  SIGMA 

DO  113  I  =  1,NS 

112  WRITE(6,511)  I 

CALL  READR  (SIGMA(D) 

IF  (SIGMA(I)  .LT.  0.)  WRITE  (6,512) 

IF  (SIGMA(I)  .LT.  0.)  GO  TO  112 

113  CONTINUE 

IF  (  CHANGE  .EQ.  NSCHAR  )  GOTO  200 
C 
C  GET  THE  NUMBER  OF  CONSTANT  WN  CURVES 

114  CONTINUE 

CALL  EXCMS  CCLRSCRN1) 

WRITE(6,513) 

CALL  READI  (NW) 

IF  (NW  .LT.  1)  GOTO  116 
C 
C  GET  THE  WN  VALUES 

WNMAX  =  WNX10XXND 

DO  116  I  =  1,NW 

115  WRITE(6,514)  I 
CALL  READR  (W(D) 

IF  (W(I)  .LT.  WN  .OR.  W(I)  .GT.  WNMAX)  WRITE  (6,515)  WN, WNMAX 
IF  (W(I)  .LT.  WN  .OR.  W(I)  . GT .  WNMAX)  GO  TO  115 

116  CONTINUE 

IF  (  CHANGE  .EQ.  NWCHAR  )  GOTO  200 
C 
C  GET  THE  NUMBER  OF  CONSTANT  ZETAXWN  CURVES 

117  CONTINUE 

CALL  EXCMS  CCLRSCRN') 

WRITE(6,516) 

CALL  READI  (NZW) 

IF  (NZW  .LT.  1)  GOTO  119 
C 
C  GET  THE  Z*WN  VALUES 

DO  119  I  =  1,NZW 

118  WRITE(6,517)  I 
CALL  READR  (ZW(D) 

IF  (ZW(I)  .LE.  0.)  WRITE  (6,518) 
IF  (ZW(I)  .LE.  0.)  GO  TO  118 


119   CONTINUE 

IF  (  CHANGE  .EQ.  ZWCHAR  )  GOTO  200 
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c 

C  GET  CONSTANT  COEFFICIENTS 

120  CONTINUE 

CALL    EXCMS    CCLRSCRN') 
WRITE(6,519) 
DO    121    I    =    NCL-1 

II    =    1-1 

WRITE(6,520)    II, I 

CALL    READR    (CJ(D) 

121  CONTINUE 

IF  (  CHANGE  .EQ.  CJCHAR  )  GOTO  200 
C 
C  GET  ALPHA  COEFFICIENTS 

122  CONTINUE 

CALL  EXCMS  CCLRSCRN') 
WRITE(6,521) 
DO  123  I  =  NC,1,-1 

II  =  1-1 

WRITE(6,522)  II, I 

CALL  READR  (AJ(D) 

123  CONTINUE 

IF  (CHANGE  .EQ.  AJCHAR  )  GOTO  200 
C 
C  GET  BETA  COEFFICIENTS 

124  CONTINUE 

CALL  EXCMS  CCLRSCRN') 
WRITE(6,523) 
DO  125  I  =  NC,1,-1 

II  =  1-1 

WRITE(6,524)  II, I 

CALL  READR  (BJ(D) 

125  CONTINUE 
C 

C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
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c 

C  REVIEW  ENTRIES 

200  CONTINUE 

CALL  EXCMS  CCLRSCRN') 

WRITE(6,525) 

CALL  READC  (REPLY) 

IF  (REPLY  .NE.  YES)  GOTO  209 
CALL  EXCMS  CCLRSCRN') 
WRITE(6,526) 

WRITE(6,527)  (LABEL(I),  1=1,  9) 
WRITE(6,528) 

WRITE(6,529)  ND, N02, NZ, NS, NW, NZW 
WRITE(6,530)  WN 
WRITE(6,531) 
IF  (NZ  .LT.  1)  GOTO  201 

WRITE(6,532)  (ZETA(M) ,M=1 , NZ) 

GOTO  202 

201  WRITE(6,533) 

202  CONTINUE 
WRITE(6,534) 

IF  (NS  .LT.  1)  GOTO  203 

WRITE(6,532)  (SIGMA(M) ,M=1 , NS) 
GOTO  204 

203  WRITE(6,533) 

204  CONTINUE 
WRITE(6,535) 

IF  (NW  .LT.  1)  GOTO  205 

WRITE(6,532)  ( W(M) ,M=1 , NW) 
GOTO  206 

205  WRITE(6,533) 

206  CONTINUE 
WRITE(6,536) 

IF  (NZW  .LT.  1)  GOTO  207 

WRITE(6,532)  (ZW(M) ,M=1 , NZW) 
GOTO  208 

207  WRITE(6,533) 

208  CONTINUE 
WRITE  (6,576) 

CALL  EXCMS(»CLRSCRN») 

WRITE(6,537) 

WRITE(6,532)  (CJ(N) , N=NC, 1 , -1 ) 

WRITE(6,538) 

WRITE(6,532)  ( AJ(N) , N=NC, 1 , -1 ) 

WRITE(6,539) 

WRITE(6,532)  ( BJ(N) , N=NC, 1 , -1) 

WRITE(6,540) 

WRITE(6,541)  XMIN,  XMAX,  YMIN,  YMAX 

209  CONTINUE 
WRITE(6,542) 

CALL  READC  (CHANGE) 

IF  (CHANGE  .NE.  YES)  GOTO  210 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

c 
c 
c 
c 
c 
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CHANGE  ROUTINE 


210 


CALL  EXCMS  ( 
WRITE(6,543) 
CALL  READC  ( 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
IF  (CHANGE  . 
GOTO  209 
CONTINUE 
WRITE(6,544) 
CALL  READC  ( 
WRITE(6,545) 
CALL  READC  ( 


•CLRSCRN') 


CHANGE) 

EQ.  TICHAR) 
WNCHAR) 
NDCHAR) 
NOCHAR) 
NZCHAR) 
NSCHAR) 
NWCHAR) 
ZWCHAR) 
CJCHAR) 
AJCHAR) 
BJCHAR) 
ENCHAR) 
PRCHAR) 
NCCHAR) 


EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 
EQ. 


TABLE) 
GRAPH) 


GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 
GOTO 


102 
103 
106 
107 
108 
111 
114 
117 
120 
122 
124 
405 
100 
210 


211  IF  (NZ)  213,213,212 

212  G  =  (10.**ND)*x(l./299.) 

213  CONTINUE 
Z  =  1.0E-60 
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c 

C  CONSTANT  ZETA  PLOTS 

IF  (NZ)  309,309,300 
300       CALL  EXCMS  ('CLRSCRN') 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  . EQ .  YES)  WRITE  (6,546) 


C 
C 


C 


DO  308  M=1,NZ 

IF  (TABLE  .EQ.  YES)  WRITE  (6,547) 
J  =  0 
R  =  0. 
WNA  =  WN 

DO  306  L=l,300 
Dl  =  0.0 
D2  =  0.0 
CI  =  0.0 
C2  =  0.0 
Bl  =  0.0 
B2  =  0.0 

DO  303  N=1,NC 
K  =  N-l 
IF  (K)  302,301,302 

301  U  =  0.0 

Ul  =  -1.0 

302  U2  =  2.0*ZETA(M)*U-U1 

Dl  =  (-1.0)XXKXCJ(N)XWNAXXKXU1+D1 

D2  =  (-1.0)*XK*CJ(N)XWNAXXKXU+D2 
CI  =  (-1.0)XXKXBJ(N)XWNAXXKXU1+C1 

C2  =  (-1.0)XXKXBJ(N)XWNAXXKXU+C2 

Bl  =  (-1.0)XXKXAJ(N)XWNAXXKXU1+B1 

B2  =  (-1.0)XXKXAJ(N)XWNAXXKXU+B2 

Ul  =  U 

U  =  U2 

303  CONTINUE 

IF  (ABS(B1XC2-B2XC1)-Z)  306,306,304 

304  J  =  J+l 

R  =  R+l. 

A(J)  =  (C1XD2-C2XD1)/(B1XC2-B2XC1) 
B(J)  =  (B2XD1-B1XD2)/(B1XC2-B2XC1) 

k 

IF  (TABLE  .NE.  YES)  GO  TO  306 

WRITE(6,548)  A(J),  B(J),  WNA,  ZETA(M) 
IF  (R/10.  -  J/10)  306,305,306 

305  CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ ,  N02) 

CALL  EXCMSCCLRSCRN') 

WRITE  (6,547) 

306  WNA  =  GXWNA 

CALL  EXCMSCCLRSCRN') 


IF  (J  .GT.  0)  GOTO  307 
WRITE(6,549) 
GOTO  308 
C 
307  IF  (GRAPH  . EQ .  YES)  CALL  PLOTD( A, B, J , . FALSE. , 

X  LABEL,  'ALPHAC, 'BETAS', 

X  MINMAX, »  Z=$',ZETA(M), 

X  XMIN,XMAX,YMIN,YMAX,GRD) 


IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1 

308  CONTINUE 

309  CONTINUE 
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c 

C  CONSTANT  SIGMA  PLOTS 

IF  (NS)  335,335,325 

325  CALL  EXCMS( 'CLRSCRN') 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  . EQ .  YES)  WRITE  (6,550) 
C 

XINC  =  (XMAX  -  XMIN)/299. 
YINC  =  (YMAX  -  YMIN)/299. 
C 

DO  334  M=1,NS 
XPT  =  XMIN 
YPT  =  YMIN 
DD  =  CJ(1) 
CC  =  BJ(1) 
BB  =  AJ(1) 
J  =  0 
R  =  0. 
C 

DO  326  N=2,NC 
K  =  N-l 
DUMMY5=SIGMA(M)XXK 

DD  =  (-1.0)XXKXCJ(N)*DUMMY5+DD 
CC  =  (-1.0)XXKXBJ(N)XDUMMY5+CC 
BB  =  (-1.0)XXKXAJ(N)XDUMMY5+BB 

326  CONTINUE 
C 

IF  (CC  .EQ.  0.  .AND.  BB  . EQ .  0.)  GOTO  334 
IF  (CC)  327,327,330 

327  DO  329  L=l,300 

J  =  J  +  l 
R  =  R+l. 
A(J)  =  XPT 

B(J)  =  (-BBXA(J)-DD)/CC 
IF  (TABLE  .NE.  YES)  GOTO  329 
WRITE(6,552)  A(J),  B(J),  SIGMA(M) 
IF  (R/10.  -  J/10)  329,328,329 

328  CALL  ROOTS  (A(J),  B(J),  A J ,  BJ,  C J ,  N02) 
CALL  EXCMS  ('CLRSCRN') 

WRITE(6,550) 

329  XPT  =  XPT  +  XINC 
GO  TO  333 

C 

330  DO  332  L=l,300 

J  =  J  +  l 
R  =  R+l. 
B(J)  =  YPT 

A(J)  =  (-CCXB(J)-DD)/BB 
IF  (TABLE  .NE.  YES)  GOTO  332 
WRITE(6,552)  A(J),  B(J),  SIGMA(M) 
IF  (R/10.  -  J/10)  332,331,332 

331  CALL  ROOTS  (A(J),  B(J),  A J ,  BJ,  CJ ,  N02) 
CALL  EXCMS  ('CLRSCRN') 

WRITE(6,550) 

332  YPT  =  YPT  +  XINC 
C 

333  CALL  EXCMS( 'CLRSCRN') 

IF  (GRAPH  .EQ.  YES)  CALL  PLOTD( A, B, J , . FALSE. , 

X  LABEL,  'ALPHAS', 'BETAS', 

X  MINMAX,'  S=$',SIGMA(M), 

X  XMIN, XMAX, YMIN, YMAX, GRD) 


IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1 

334  CONTINUE 

335  CONTINUE 
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C  CONSTANT  ZETA-OMEGA  PLOTS 

IF  (NZW)  359,359,350 
350       CALL  EXCMS('CLRSCRN') 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  . EQ .  YES)  WRITE  (6,553) 
C 


XWN 

=  WN 

DO  358  M= 

=  1, 

.NZW 

IF  (TABLE  .EQ. 

YES) 

J  =  0 

R  =  0. 

AZETA  ' 

■    ] 

L./300, 

DO  356 

L  = 

=1,300 

XWN  i 

=  ZW(M)/AZETA 

Dl  = 

0, 

,0 

D2  = 

0 

.0 

CI  = 

0 

,0 

C2  = 

0 

,0 

Bl  = 

0 

.0 

B2  = 

0 

,0 

YES)  WRITE  (6,554) 


C 


DO  353  N=1,NC 
K  =  N-l 
IF  (K)  352,351,352 

351  Ql  =  0.0 

Q  =  -1.0/XWN**2 

352  D2  =  CJ(N)*Q1+D2 
C2  =  BJ(N)*Q1+C2 
B2  =  AJ(N)*Q1+B2 
Dl  =  CJ(N)*Q+D1 
CI  =  BJ(N)*Q+C1 
Bl  =  AJ(N)*Q+B1 

Q2  =  -2.0XZW(M)*Q1-XWNXX2*Q 
Q  =  Ql 

353  Ql  =  Q2 

IF  (ABS(B1XC2-B2*C1)-Z)  356,356,354 

354  J  =  J+l 
R  =  R+l. 

A(J)  =  (C1*D2-C2XD1)/(B1XC2-B2*C1) 

B(J)  =  (B2*D1-B1*D2)/(B1XC2-B2*C1) 
C 

IF  (TABLE  .NE.  YES)  GO  TO  356 
WRITE(6,552)  A(J),  B(J),  ZW(M) 
IF  (R/10.  -  J/10)  356,355,356 

355  CALL  ROOTS  (A(J),  B(J),  AJ,  BJ,  CJ ,  N02) 
CALL  EXCMS('CLRSCRN') 

WRITE  (6,554) 

356  AZETA  =  AZETA+( 1 ./300 . ) 
C 

CALL  EXCMS('CLRSCRN') 
C 

IF  (J  .GT.  0)  GOTO  357 

WRITE(6,549) 

GOTO  358 
C 

357  IF  (GRAPH  . EQ .  YES)  CALL  PLOTD( A, B, J , . FALSE. , 

X  LABEL,  'ALPHAS', 'BETAS', 

X  MINMAX, 'ZW=$',ZW(M), 

X  XMIN,XMAX,YMIN,YMAX,GRD) 


IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1 

358  CONTINUE 

359  CONTINUE 
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c 

C  CONSTANT  OMEGA  PLOTS 

375  IF  (NW)  385,385,376 

376  CALL  EXCMS('CLRSCRN') 

IF  (TABLE  .EQ.  YES  .OR.  GRAPH  . EQ .  YES)  WRITE  (6,555)  . 
C 
C 

DO  384  M=1,NW 

IF  (TABLE  .EQ.  YES)  WRITE  (6,547) 
J  =  0 
R  =  0. 
AZETA  =0.0 
C 

DO  382  L=l,300 
Dl  =  0.0 
D2  =  0.0 
CI  =  0.0 
C2  =  0.0 
Bl  =  0.0 
B2  =  0.0 
C 

DO  379  N=1,NC 
K  =  N-l 
IF  (K)  378,377,378 

377  U  =  0.0 

Ul  =  -1.0 

378  U2  =  2.0*AZETA*U-U1 

Dl  =  (-1.0)XXK*CJ(N)XW(M)X*K*U1+D1 
D2  =  (-1.0)XXKXCJ(N)XW(M)XXKXU+D2 
CI  =  (-1.0)XXKXBJ(N)XW(M)X*KXU1+C1 
C2  =  (-1.0)XXKXBJ(N)XW(M)XXKXU+C2 
Bl  =  (-1.0)XXKXAJ(N)XW(M)XXKXU1+B1 
B2  =  (-1.0)XXKXAJ(N)XW(M)XXKXU+B2 
Ul  =  U 
U  =  U2 

379  CONTINUE 

IF  (ABS(B1XC2-B2XC1)-Z)  382,382,380 

380  J  =  J+l 

R  =  R+l. 

A(J)  =  (C1XD2-C2XD1)/(B1XC2-B2XC1) 
B(J)  =  (B2XD1-B1XD2)/(B1XC2-B2XC1) 
C 

IF  (TABLE  .NE.  YES)  GO  TO  382 

WRITE(6,548)  A(J),  B(J),  W(M),  AZETA 

IF  (R/10.  -  J/10)  382,381,382 

381  CALL  ROOTS  (A(J),  B(J),  A J ,  BJ,  CJ ,  N02) 
CALL  EXCMS( 'CLRSCRN') 

WRITE  (6,547) 

382  AZETA  =  AZETA+( 1 ./299  .  ) 
CALL  EXCMS('CLRSCRN') 

C 

IF  (J  .GT.  0)  GOTO  383 

WRITE(6,549) 

GOTO  384 
C 

383  IF  (GRAPH  . EQ .  YES)  CALL  PLOTD(A, B, J , . FALSE. , 

X  LABEL,  'ALPHAS', 'BETAS', 

X  MINMAX,'  W=$',W(M), 

X  XMIN,XMAX,YMIN,YMAX,GRD) 

C 

384  CONTINUE 
C 

385  CONTINUE 
C 

C 
C 
C 
C 
C 
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c 


IF  (GRAPH  .EQ.  YES)  CALL  PLOTD( 0 . , 0 . , 0 . , .TRUE. , 
X  •  *»,•  $»,«  $«, 

X  0.,'  $',-9.7531, 

X  0.,0.,0.,0.,GRD) 

IF  (GRAPH  .EQ.  YES)  GRD  =  GRD+1 . 
GRD  =  0. 

400  CALL  EXCMSCCLRSCRN') 

IF  (OPT  .EQ.  YES)  CALL  DONEPL 

IF  (OPT  .EQ.  YES)  STOP 

WRITE(6,557) 

WRITE(6,556) 

WRITE(6,557) 

WRITE(6,558) 

WRITE(6,559) 

WRITE(6,560) 

WRITE(6,561) 

WRITE(6,562) 

WRITE(6,563) 

WRITE(6,557) 

CALL  READI  (IANS) 

IF  (IANS  .GT.  6  .OR.  IANS  .LT.  1)  GOTO  400 
GO  TO  (100,  404,  401,  402,  403,  405)  IANS 
C 
C  ROOT  FINDER  OPTION 

401  CALL  EXCMSCCLRSCRN*  ) 
WRITE(6,564) 
WRITE(6,565) 

CALL  READR  (ALPHA) 

WRITE(6,566) 

CALL  READR  (BETA) 

CALL  ROOTS  (ALPHA,  BETA,  AJ,  B J ,  CJ ,  N02) 

GO  TO  400 
C 
C  SAVE  OPTION 

402  CALL  SAVIT 
GO  TO  400 

C 

C  CREATE  DISSPLA  METAFILE  OPTION 

403  CALL  EXCMSCCLRSCRN') 
WRITE(6,567) 
WRITE(6,568) 
WRITE(6,569) 

CALL  READC  (OPT) 

IF  (OPT  .NE.  YES)  GOTO  400 

CALL  EXCMSCCLRSCRN' ) 

WRITE(6,570) 

CALL  DONEPL 
C 
C  COMPRS  =  SUBROUTINE  TO  LET  DONEPL  FINISH 

CALL  META 

GO  TO  211 
C 
C  SAME  PROBLEM  OPTION 

404  WRITE(6,571) 

CALL  READI  (MINMAX) 

IF  (MINMAX  .EQ.  1)  GO  TO  200 

CALL  EXCMSCCLRSCRN') 

WRITE(6,572) 

CALL  READR  (XMIN) 

WRITE(6,573) 

CALL  READR  (XMAX) 

WRITE(6,574) 

CALL  READR  (YMIN) 

WRITE(6,575) 

CALL  READR  (YMAX) 
GO  TO  200 

405  CONTINUE 
RETURN 


102 


500 


501 
502 
503 
504 
505 
506 
507 
508 
509 
510 

511 
512 
513 
514 
515 


516 
517 
518 
519 
520 
521 
522 
523 
524 
525 
526 
527 
528 
529 
530 
531 
532 
533 
534 
535 
536 
537 
538 
539 
540 
541 
542 
543 


FORMAT 

X 

X 

X 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

X  (NS)' 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

X 

X 

FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 

X5X,'TI 


FORMATS  STATEMENTS 
(5X,'THIS  IS  THE  INTERACTIVE  PARAMETER  PLANE  PROGRAM...',/, 

5X,'THE  USER  WILL  BE  PROMPTED  FOR  VARIOUS  INPUTS.',////, 

5X,*WILL  YOU  BE  ENTERING  DATA  FROM  A  CONSOLE  OR  DATAFILE?', 

/  15X  '"D"   OR   "C" ' ) 
(///, ''ENTER  TITLE  TO  APPEAR  FOR  THIS  FAMILY  OF  CURVES') 
C/, IX, 'WHAT  IS  THE  STARTING  VALUE  OF  OMEGAN  (WN>0.0)?') 
(/,1X,'WN  MUST  BE  GREATER  THAN  ZERO  -  TRY  AGAIN') 
(/,1X,'H0W  MANY  DECADES  PAST  WN  ARE  DESIRED?  (ND)') 
(/,'  WHAT  IS  THE  ORDER  OF  THE  CHARACTERISTIC  EQUATION?( NO) ' ) 
V, IX, 'HOW  MANY  CONSTANT  ZETA  CURVES  ARE  DESIRED?  (NZ)') 
(/,'  ENTER  THE  VALUES  OF  ZETA  TO  BE  USED  IN  COMPUTATION...') 
(5X,  'ZETAC',12,  ')  =  ?») 

(5X,»ZETA  MUST  LIE  BETWEEN  0  AND  1,  INCLUSIVE  -  TRY  AGAIN') 
(/,'  HOW  MANY  CONSTANT  SIGMA  CREAL  ROOT)  CURVES  ARE  DESIRED? 
) 

(5X,  'SIGMACI2,  ')  =  ?') 

(5X, 'NEGATIVE  SIGMA  MEANS  POSITIVE  REAL  ROOT  -  TRY  ANOTHER') 
(/,1X,'H0W  MANY  CONSTANT  WN  CURVES  ARE  DESIRED?  (NW)«) 
(5X, »W(',I2, ')  =  ?') 

(5X,'WN  NOT  WITHIN  PLOTTABLE  RANGE', 
/,5X,'Y0UR  USABLE  RANGE  IS' 
/,10X,F10.2, '   TO  ',F10.2) 

(/,1X,'H0W  MANY  CONSTANT  Z*WN  CURVES  ARE  DESIRED?  (NZW)') 
(5X,  'ZWCI2,  »)  =  ?») 

(5X, 'NON-POSITIVE  Z-WN  MEANS  POSITIVE  ROOT  -  TRY  ANOTHER') 
(/, IX, 'ENTER  THE  CONSTANT  COEFFICIENTS...') 
C5X,  ' SXX',12,' CJCI2,')  =  ?') 


(/, IX, 'ENTER  THE  ALPHA  COEFFICIENTS...') 

C5X, ' SXXM2,  ■ AJCI2,')  =  ?») 

(/, IX, 'ENTER  THE  BETA  COEFFICIENTS...') 

( 5X  i s*x'  12  ' BJ('  12  ')  =  7') 

(/,'  WANT  TO' REVIEW  YOUR  ENTRIES  BEFORE  RUNNING?  CY/N)') 

(/,10X, '  GRAPH  TITLE') 

(1X,9A4) 

(/,8X,2HND,8X,2HN0,8X,2HNZ,8X,2HNS,8X,2HNW,7X,3HNZW) 

(6110) 

(/,10X, 'INITIAL  VALUE  OF  OMEGA  =  »,F10.5) 

(/,10X, '  ZETA  ») 

(8E10.3) 

(IX,  ' NO  VALUE ') 

(/,10X, 'SIGMA  ') 

(/,10X,'  W     ») 

(/,10X,'ZW     ') 

(/,10X, 'CONSTANT  COEFFICIENTS  IN  DECENDING  ORDER') 

(/,10X, 'ALPHA     COEFFICIENTS  IN  DECENDING  ORDER') 

(/,10X,'BETA      COEFFICIENTS  IN  DECENDING  ORDER') 

(/,10X,'XMIN    XMAX    YMIN    YMAX') 

(1X,4E10.3) 
(/,'  WANT  TO  MAKE  ANY  CHANGES?  (Y/N)') 
(/,»  WHAT  VARIABLE/AREA  DO  YOU  WISH  TO  CHANGE?',//, 
TLE TI      OMEGA  START.  .WN      #  DECADES ND',/ 
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c 
c 
c 
c 
c 
c 
c 
c 
c 


544 
545 
546 
547 
548 
549 

550 
551 
552 
553 
554 
555 
556 
557 
558 
559 
56  0 
561 
562 
563 
564 

565 
566 
567 


568 

569 
57  0 

571 

572 
57  3 
574 
575 
576 


C 
C 
C 
C 
C 
C 

c 
c 
c 
c 


.BJ',/ 
.NC',///, 


X5X, 'ORDER NO      CONST  ZETA ... NZ      CONST  SIGMA .. NS ' ,/ 

X5X, 'CONST  WN NW      CONST  Z-WN . . .ZW* ,// 

X5X, 'CONST  TERMS.. CJ      ALPHA  TERMS.. AJ      BETA  TERMS 

X5X,'END E      NEW  PROBLEM.  .PR      NO  CHANGE.  . 

XIX, 'ENTER  TWO  DIGIT  CODE...') 

FORMATC/,'  DO  YOU  WANT  TABULATED  DATA  ON  THE  SCREEN?  (Y/N)») 

FORMATC/,'  DO  YOU  WANT  THE  GRAPHS  ON  THE  TERMINAL?  CY/N)') 

FORMAT  C1H1,10X, 'CONSTANT  ZETA  CURVES') 

FORMATC/, 10X, 'ALPHA  BETA  OMEGA  ZETA') 

FORMAT  (4E16.5) 

FORMATC/,'  DUE  TO  PLOT  RESTRICTIONS,  A  COMPLETE  GRAPH  CANNOT  BE  OU 
XTPUT. •) 

FORMAT  (1H1,10X, 'CONSTANT  SIGMA  CURVES') 

FORMAT  (/,10X, 'ALPHA  BETA  SIGMA') 

FORMAT  (3E16.5) 

FORMAT  (1H1,10X, 'CONSTANT  ZETA-OMEGA  CURVES') 

FORMAT  (/,10X, 'ALPHA  BETA         ZETA-OMEGA') 

FORMAT  (1H1,10X, 'CONSTANT  OMEGA  CURVES') 

F0RMAT(5X,'|  OPTION  NO. I  OPTION 

F0RMATC5X, 

F0RMAT(5X, 

F0RMATC5X, 

F0RMATC5X, 

F0RMAT(5X, 

F0RMATC5X, 

F0RMATC5X, 

F0RMAT(5X, 
X        5X 

FORMATC/// 


FORMATC/, 5X, 'ENTER  THE  BETA  VALUE 


FORMATC 5X, 
X  5X, 
X  5X, 
X  5X, 
X        5X 

F0RMATC5X! 
X  5X, 
X        5X 

FORMATC/// 
X        15X 


NEW  PROBLEM 

SAME  PROBLEM 

ROOT  FINDER 

SAVE  DATA 

SAVE  GRAPH  IN  DISSPLA  METAFILE 

RETURN  TO  MAIN  MENU 
ENTER  AN  ALPHA-BETA  PAIR,  AND  THE  ROOTS  OF  YOUR 
SYSTEMS  CHARACTERISTIC  EQUATION  WILL  BE  RETURNED 
5X, 'ENTER  THE  ALPHA  VALUE 


YOU  NOW  HAVE  THE  OPTION  OF  STORING  THE  LAST  SET  OF   »/ 

CURVES  IN  A  DISSPLA  METAFILE.   THIS  ALLOWS  RETRIEVAL'/ 

OF  DATA  AT  A  LATER  TIME  FOR  ROUTING  TO  ANY  OF        •/ 

SEVERAL  OUTPUT  DEVICES  CTEK618,  3800  LASER  PRINTER,  '/ 

VERSATEC  PLOTTER,  ETC.) 

IF  YOU  CHOOSE  THIS  OPTION,  THE  PROGRAM  MUST  BE 

TERMINATED  -  THIS  CANNOT  BE  AVOIDED  WITHOUT 

CATASTROPHIC  RESULTS. 

5X,'D0  YOU  WISH  TO  USE  THIS  OPTION? 

myn   QR   "N" 


FORMATC/////, 5X, 'IF  YOU  WISH  GRAPHIC  OUTPUT,  TYPE: 

X        15X, '"DISSPOP" 

X        5X,'AND  FOLLOW  THE  INSTRUCTIONS  .  .  . 
FORMAT  C///,'  AUTOSCALE  OR  USER-DEFINED  LIMITS  FOR  CURVES?', 

X/,'  l=AUTOSCALE;  0=USER-DEFINED' ) 
FORMAT  C//,»  INPUT  MINIMUM  VALUE  FOR  X  CX-MIN)*) 
FORMAT  C/,'  INPUT  MAXIMUM  VALUE  FOR  X  CX-MAX)') 
FORMAT  C/,'  INPUT  MINIMUM  VALUE  FOR  Y  CY-MIN)') 
FORMAT  C/,'  INPUT  MAXIMUM  VALUE  FOR  Y  CY-MAX)') 
FORMAT  C////////////////) 
END 


»// 
'/ 
»/ 
'/// 
'/ 
'/ 
'/ 


»/ 
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SUBROUTINE  PLOTD  —  GRAPHS  WILL  BE  PRODUCED  ON  UPRIGHT  11  X  14 
PAGE  WITH  9  INCH  AXES.   IF  USER  SELECTS  'AUTOSCALE'  FEATURE, 
SUBROUTINE  PLD010  (INTERNAL  TO  PLOTD)  FINDS  MIN  AND  MAX  FOR  EACH 
AXIS  AND  SCALES  ACCORDINGLY.   FORMAT: 

CALL  PLOTDCXDATA,  YDATA,  NNPTS,  EJECT,  LABEL,  XLABEL,  YLABEL, 
X  MINMAX,  CRVTTL,  CRVNUM,  XMIN,  XMAX,  YMIN,  YMAX) 

WHERE: 

XDATA      IS  A  REAL**  ARRAY  DIMENSIONED  AT  LEAST  I NNPTS I 
CONTAINING  THE  X  ORDINATE  VALUES, 

YDATA      IS  A  REAL**  ARRAY  DIMENSIONED  AT  LEAST  I NNPTS I 
CONTAINING  THE  Y  ORDINATE  VALUES, 

NNPTS      IS  AN   INTEGER**  SCALAR  DESIGNATING  THE   NUMBER  OF 
POINTS  TO  BE  PLOTTED.   THE  NUMBER  OF  POINTS  IS 
ABS(NNPTS).   NNPTS<0  MEANS  PLOT  POINTS  ONLY. 

EJECT  IS  A  LOGICAL**  VARIABLE  OR  CONSTANT  INDICATING 
WHETHER  A  PAGE  EJECT  IS  REQUIRED  FOLLOWING  THE 
CURRENT  CURVE.  THIS  ALLOWS  MULTIPLE  CURVES  ON  ONE 
SET  OF  EXES.  PAGE  EJECT  WILL  OCCUR  FOR  NEXT  GRAPH 
AFTER  EJECT  HAS  BEEN  SET  TO  .TRUE. 

LABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY 
CONTAINING  THE  INTENDED  LABEL  FOR  THE  GRAPH.  THE 
MAXIMUM  ALLOWABLE  LENGTH  (INCLUDING  '$»  CHARACTER) 
IS  32  CHARACTERS. 

XLABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY 
CONTAINING  THE  INTENDED  LABEL  OF  THE  X-AXIS  OF  THE 
GRAPH.   IN  THIS  PROGRAM,  XLABEL  IS  ALWAYS  'ALPHA'. 

YLABEL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY 
CONTAINING  THE  INTENDED  LABEL  OF  THE  Y-AXIS  OF  THE 
GRAPH.   IN  THIS  PROGRAM,  YLABEL  IS  ALWAYS  'BETA'. 

MINMAX  IS  A  PARAMETER  THAT  DETERMINES  WHETHER  THE  MINIMUM 
AND  MAXIMUM  VALUES  FOR  THE  AXES  ARE  TO  BE  ASSIGNED 
BY  THE  USER,  OR  WHETHER  THEY  WILL  BE  ' AUTOSCALED' . 

CRVTTL  IS  A  QUOTED  LITERAL  OR  HOLLERITH  STRING  OR  ARRAY 
AND  TERMINATED  BY  A  ■$•  CHARACTER  SPECIFYING  THE 
INTENDED  NAME  WHICH  LABELS  AN  INDIVIDUAL  CURVE. 

CRVNUM  IS  A  REAL  VARIABLE  OR  CONSTANT  THAT  SPECIFIES  THE 
VALUE  TO  BE  CONCATENATED  ONTO  THE  END  OF  'CRVTTL'. 
FOR  EXAMPLE,  IF  THIS  CURVE  REPRESENTS  'ZETA  =  0.5' 
THEN  CRVTTL  =  'Z=$',  WHILE  CRVNUM  =  0.5. 


SUBROUTINE  PLOTDCXDATA,  YDATA,  NNPTS,  EJECT,  LABEL,  XLABEL, 

X  YLABEL,  MINMAX,  CRVTTL,  CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX,GRD) 

REAL**  XDATA(l),  YDATA(l) 

REAL**  XMIN,  XMAX,  YMIN,  YMAX 

REAL**  CRVTTL, CRVNUM 

INTEGER**  NPTS, NNPTS 

LOGICAL**  EJECT 

LOGICAL*l  LABEL(l) 

LOGICAL*l  XLABEL(l) 

LOGICAL*!  YLABEL(l) 


LOGICAL** 


INIT  /.FALSE./ 
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c 
c 
c 

C      SET  THE  NUMBER  OF  POINTS 
C 

NPTS  =  IABS(NNPTS) 
C 

C         IF  THE  ROUTINE  HAS  BEEN  NOT  BEEN  INITIALIZED  (INIT  =  .FALSE.) 
C        THEN  INITIALIZE  IT. 
C 

IF  (.NOT.  INIT)  CALL  PLDOOKXDATA,  YDATA,  NPTS,  LABEL,  XLABEL, 
X  YLABEL,  MINMAX,  CRVTTL , CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX,GRD) 

C 

C      FRAME  THE  PLOT  AND  REORDER  SYMBOLS 
C 

IF  (.NOT.  INIT)  CALL  FRAME 
C 

C         INDICATE  INITIALIZATION  IN  CASE  THIS  ROUTINE  IS  RE-ENTERED 
C        WITH  MULTIPLE  CURVES. 
C 

INIT  =  .TRUE. 
C 

C        CALCULATE  THE  NUMBER  OF  POINTS  TO  GIVE  ABOUT  5  MARKERS  PER 
C         LINE 
C 

C      TENTATIVELY  SET  POINTS  ONLY,  THEN  CHECK 
NMARK  =  -1 

IF(NNPTS.LE.0)GO  TO  10 
c 
C      CURVE  WANTED,  SET  MARKERS 

IF  (MOD(NPTS,  4)  .EQ.  1)  NMARK  =  NPTS  /  4 
IF  (MOD(NPTS,  4)  .NE.  1)  NMARK  =  NPTS  /  3 
IF  (NMARK  .EQ.  0)  NMARK  =  1 
10   CONTINUE 
C 
C 

C         DRAW  THE  CURVE 
C 

XXMAX  =  -1.0E75 
YYMAX  =  -1.0E75 
C 

DO  20  1=1, NPTS 

IF(XDATAd)  .GT.  XMAX  .OR.  XDATA(I)  .LT.  XMIN  .OR. 
X      YDATA(I)  .GT.  YMAX  .OR.  YDATA(I)  .LT.  YMIN)  GO  TO  20 
IF(XXMAX  .LT.  XDATA(D)  J  =  I 
IF(XXMAX  .LT.  XDATA(D)  XXMAX  =  XDATA(I) 
IF(YYMAX  .LT.  YDATA(D)  K  =  I 
IF(YYMAX  .LT.  YDATA(D)  YYMAX  =  YDATA(I) 
20   CONTINUE 
C 

IFUYMAX-YYMAX)  -  (XMAX-XXMAX) )  40,30,30 
30      L  =  J 

GO  TO  50 
40      L  =  K 
c 

50   CALL  GRACE(0.) 
CALL  DOT 

IF  (GRD)  60,60,70 
60   CALL  GRID(1,1) 
70   CALL  RESET( 'DASH') 

IF  (CRVNUM  .EQ.  -9.7531)  GO  TO  80 
CALL  HEIGHT(0.125) 

CALL  RLMESS(CRVTTL,3,XDATA(L),YDATA(D) 
CALL  RLREAL(CRVNUM,2, 'ABUT', 'ABUT') 
CALL  THKCRV(0.015) 
80   CALL  CURVE(XDATA,  YDATA,  NPTS,  NMARK) 
CALL  RESET('THKCRV') 
C 
C        IF  THIS  IS  NOT  THE  LAST  (OR  ONLY)  CURVE  ON  THIS  GRAPH,  THEN 
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C         EXIT.   OTHERWISE  CLOSE  THE  PLOT  AND  TURN  OFF  INITIALIZATION 

C         FLAG. 

C 

IF  (  .NOT.  EJECT)  RETURN 
C 

C         END  OF  THIS  PLOT 
C 

CALL  ENDPL(O) 

INIT  =  .FALSE. 

RETURN 

END 


C 
C 

SUBROUTINE  PLDOOKXDATA,  YDATA,  NPTS,  LABEL,  XLABEL, 
X  YLABEL,  MINMAX,  CRVTTL,  CRVNUM, 

X  XMIN,  XMAX,  YMIN,  YMAX,  GRD) 

C 

C        THIS  SUBROUTINE  DOES  THE  INITIALIZING. 
C 

REAL*4        XDATA(NPTS) 

REAL**        YDATA(NPTS) 

INTEGER**     NPTS 

LOGICALXl     LABEL(l) 

LOGICALXl     XLABEL(l) 

LOGICALXl     YLABEL(l) 
C 
C 

REALX4        XMIN 

REALX4        XMAX 

REALX4        YMIN 

REALX4        YMAX 
C 

C         INITIALIZE  DISSPLA 
C 

CALL  PLD009 

CALL  HEADIN(LABEL,  100,  2.,  1) 

CALL  XNAMECXLABEL,  100) 

CALL  YNAME(YLABEL,  100) 
C 
C 

C         EXTRACT  MINIMA  AND  MAXIMA 
C 

IF  (MINMAX  .NE.  1)  GO  TO  90 

CALL  PLD010(XDATA,  NPTS,  XMIN,  XMAX) 

CALL  PLD010(YDATA,  NPTS,  YMIN,  YMAX) 
C 

C        CALL  THE  LINEAR-LINEAR  INITIALIZING  ROUTINE 
C 

90   CALL  PLD011  (XMIN, XMAX, YMIN, YMAX) 

RETURN 

END 
C 
C 

SUBROUTINE  PLD009 
C 

C         THIS  SUBROUTINE  ESTABLISHES  THE  PARAMETERS  FOR  DISSPLA. 
C 

C         NOTE  THAT  IT  IS  THE  USER'S  RESPONSIBILITY  TO  NOMINATE  THE 
C        GRAPHIC  DEVICE. 
C 
C 

CALL  NOCHEK 

CALL  NOBRDR 

CALL  PAGE(14.,14.) 

CALL  PHYS0R(2., .75) 
C 

C        GO  TO  LEVEL  2. 
C 

CALL  AREA2D(9.,9.) 
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C         LETTERING  IS  DUPLX  WITH  UPPER  CASE  ONLY. 
C 

CALL  DUPLX 

CALL  BASALF( 'STANDARD') 
C 

C         INTEGER  (OR  ROUNDED)  AXES  WITH  Y  AXIS  LABELLING  AT  0  DEGREES. 
C 

CALL  INTAXS 

CALL  YAXANGCO.) 

RETURN 

END 
C 
C 

SUBROUTINE  PLD010CV,  N,  MIN,  MAX) 
C 

C        THIS  SUBROUTINE  SCANS  A  VECTOR  FOR  MAXIMUM  AND  MINIMUM 
C 

C         INPUT  PARAMETERS? 
C 

C  V  DATA  VECTOR  (REAL) 

C  N  NUMBER  OF  POINTS  IN  VECTOR  V  (INTEGER) 

C 

C         OUTPUT  PARAMETERS: 
C 

C  MIN        VECTOR  ORIGIN  (REAL) 

C  MAX        VECTOR  MAXIMUM  (REAL) 

C 

REAL**        V(N) 

INTEGER**     N 

REAL**        MIN 

REAL**        MAX 
C 

C         INITIALIZE  THE  MAXIMA  AND  MINIMA 
C 

MIN  =  1.0E75 

MAX  =  -1.0E75 
C 

C         FIND  MAXIMUM  AND  MINIMUM  OF  VECTOR  V 
C 

DO  100  I  =  1,  N 

IF  (MIN  .GT.  V(D)  MIN  =  V(I) 

IF  (MAX  .LT.  V(D)  MAX  =  V(I) 
100   CONTINUE 

RETURN 

END 
C 
C 

SUBROUTINE  PLD011  (XMIN, XMAX, YMIN, YMAX) 
C 

C        THIS  SUBROUTINE  SETS  UP  DISSPLA  FOR  A  LINEAR-LINEAR  AXIS  PLOT 
C 

REAL**        XMIN 

REAL**        XMAX 

REAL**        YMIN 

REAL**       YMAX 
C 

C        A  SIMPLE  CALL  TO  GRAF  WILL  DO  IT... 
C 

CALL  HEIGHT(0.175) 

CALL  GRAF(XMIN,  'SCALE',  XMAX,  YMIN,  'SCALE',  YMAX) 

RETURN 

END 
C 
C 
C 
C 
C 
c 
c 
c 
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c 

C  ================================================================= 

C   SUBROUTINE  ROOTS  ~  CALCULATES  ROOTS  OF  THE  NO  ORDER  EQUATION 

C   FOR  EVERY  TENTH  VALUE  ALPHA/BETA  PAIR. 

C  ================================================================= 

c 
c 

SUBROUTINE  ROOTS  (ALPHA,  BETA,  AJ,  B J ,  CJ ,  N02) 
INTEGER**  N02,  NN,  NNN,  NNNN 

REALX4  ALPHA,  BETA,  AJ(IOO),  BJ(IOO),  CJ(IOO) 
REAL*8  COEF(IOO),  ROOTMY(IOO) 
WRITE(6,<+0) 
NN  =  N02+1 
NNNN  =  N02 
DO  10  1=1, NN 
NNN  =  NN+l-I 
10   COEF(NNN)  =  CJ(I)  +  (ALPHA  x  AJ(D)  +  (BETA  x  BJ(D) 
CALL  ZPOLR  (C0EF,N02,R00TMY,IER) 
WRITE  (6,50) 
DO  20  1=1, NNNN 

II  =  2*1 

III  =  II-l 

20   WRITE  (6,60)  ROOTMY(III),  ROOTMY(II) 

WRITE  (6,30) 
30   FORMAT(////////////////////) 
40   FORMAT(/20X, »  ROOTS  FOR  ABOVE  ALPHA,  BETA*) 
50   F0RMAT(21X,'  REAL  PART        IMAG  PART') 
60   F0RMAT(21X,E10.4,6X,E10.4) 
RETURN 
END 
C 
C 
C  ================================================================= 

C   SUBROUTINE  SAVIT  ~  SAVES  DATA  IN  FN  FT  FM  =  INAME  DATA  Al, 

C   WHERE  INAME  IS  THE  USER'S  CHOICE. 

C  ================================================================= 

c 
c 

SUBROUTINE  SAVIT 

COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  C J ,  AJ,  BJ, 
X  NZ,  ZETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

REALX4  WN,  CJ(IOO),  AJ(IOO),  BJ(IOO),  XMIN,  XMAX,  YMIN,  YMAX 

REALX4  ZETA(IOO),  SIGMA(IOO),  W(100),  ZW(IOO) 

INTEGER  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 

CHARACTERX4  LABEL(9),  INAME(2) 

WRITE(6,10) 

READ(5,20)  (INAME(I),  1=1,2) 

CALL  FRTCMSCFILEDEF  ','02       ','DISK     ',  INAME,  ■  DATA     ' 

WRITE(2,30)  (LABEL(I),  1=1,  9) 

WRITE(2,x)  WN 

WRITE(2,X)  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 

WRITE(2,X)  (CJ(J),  J=l,  NC) 

WRITE(2,X)  (AJ(J),  J=l,  NC) 

WRITE(2,x)  (BJ(J),  J=l,  NC) 

WRITE(2,x)  (ZETA(M),  M=l,  NZ) 

WRITE(2,*)  (SIGMA(M),  M=l,  NS) 

WRITE(2,x)  (W(M),  M=l,  NW) 

WRITE(2,x)  (ZW(M),  M=l,  NZW) 

WRITE(2,x)  XMIN,  XMAX,  YMIN,  YMAX 
10   F0RMAT(5X, 'UNDER  WHAT  NAME  DO  YOU  WANT  TO  SAVE  THE  DATA?',/, 

X        5X,»(8  CHARACTERS  MAX)') 
20   F0RMAT(2A4) 
30   F0RMAT(9A4) 

END  FILE  02 

REWIND  02 

RETURN 

END 
C 

c 
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c 

c  ================================================================== 

C   SUBROUTINE  GETIT  —  RETRIEVES  DATA  FROM  FN  FT  FM  =  INAME  DATA  Al, 

C   WHERE  INAME  IS  THE  USER'S  CHOICE. 

C  ================================================================== 

c 
c 

SUBROUTINE  GETIT 

COMMON  /SAVE/  LABEL,  WN,  ND,  N02,  NC,  C J ,  AJ,  B  J , 
X  NZ,  ZETA,  NS,  SIGMA,  NW,  W,  NZW,  ZW, 

X  XMIN,  XMAX,  YMIN,  YMAX 

REAL**  WN,  CJ(IOO),  AJ(IOO),  BJ(IOO),  XMIN,  XMAX,  YMIN,  YMAX 
REAL*4  ZETA(IOO),  SIGMAC100),  W(100),  ZWC100) 
INTEGER  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 
CHARACTER**  LABELC9) 
CHARACTERX8  INAME 
CHARACTERX21  NAME 
10   WRITE(6,40) 

READ(5,50)  INAME 

NAME  =  'STATE  '//INAME//'  DATA  X» 

CALL  EXCMS(NAME,RC) 

IF  (RC  .EQ.  0)  GOTO  20 
WRITE(6,30) 
GOTO  10 
20   CALL  FRTCMSCFILEDEF  ','02       ','DISK     ',  INAME,  '  DATA     ') 
READ(2,60)  (LABEL(I),  1=1,  9) 
READ(2,X)  WN 
READ(2,X)  ND,  N02,  NC,  NZ,  NS,  NW,  NZW 


READ(2,X)  (CJ(J),  J=l,  NC) 

READ(2,X)  (AJ(J),  J=l,  NC) 

READ(2,X)  (BJ(J),  J=l,  NC) 

READ(2,X)  (ZETA(M),  M=l,  NZ) 

READ(2,X)  (SIGMA(M),  M=l,  NS) 

READ(2,X)  (W(M),  M=l,  NW) 

READ(2,X)  (ZW(M),  M=l,  NZW) 

READ(2,X)  XMIN,  XMAX,  YMIN,  YMAX 

30 

F0RMATC5X, 'DATA  FILE  NOT  FOUND  - 

40 

F0RMATC5X, 'UNDER  WHAT  NAME  IS  THE 
X        5X,'(S  CHARACTERS  MAX)') 

50 

F0RMATUA8) 

60 

F0RMAK9A4) 
END  FILE  02 
REWIND  02 
RETURN 

C 
C 

END 

TRY  ANOTHER') 
DATAFILE  SAVED?        ',/, 


C  ============================================================= 

C  SUBROUTINE  ASTER  PLACES  A  $  AT  THE  END  OF  A  CHARACTER  STRING 
C  ============================================================= 

c 
c 

SUBROUTINE  ASTERCLLINES, LINE) 

CHARACTERX4  DOLLAR/'$    '/,  LLINESC8),  LINEC9) 

DO  10  1=1,8 

LINE(I)=LLINES(I) 
10     CONTINUE 

LINE(9)=DOLLAR 

RETURN 

END 
C 
C 
C 
C 
C 
C 
C 
C 

c 
c 
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c 

c  ==================================================================== 

C   SUBROUTINE  READR  ~  INTERACTIVELY  READS  A  REAL  NUMBER  REPLY. 

C   IF  THE  USER  INADVERTENTLY  ENTERS   A  NULL  STRING 

C   A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 

C  ==================================================================== 

c 
c 
c 

SUBROUTINE  READR  (ANSR) 

REALX4  ANSR 

INTEGER  COUNT 

COUNT=0 
10     CONTINUE 

COUNT=COUNT+l 

IF  (COUNT. LT. 3)  GO  TO  20 

WRITE  (5,60) 

GO  TO  40 
20     CONTINUE 

READ  (5,*,END=30,ERR=30)  ANSR 

RETURN 
30     REWIND  5 

WRITE  (5,50) 

GO  TO  10 
40     CONTINUE 

STOP 
50     FORMAT  (IX,'  WARNING:   NULL  STRINGS  ARE  NOT  ALLOWED,  ENTER  A  NUMER 

XICAL  VALUE. •) 
60     FORMAT  (///,5X,*  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 

END 
C 
C 
C  ==================================================================== 

C   SUBROUTINE  READI  —  INTERACTIVELY  READS  AN  INTEGER  REPLY. 

C   IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING  OR  NEGATIVE  VALUE   = 

C   A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 

C  ==================================================================== 

c 
c 

SUBROUTINE  READI  (IANS) 

INTEGER  COUNT, IANS 

COUNT=0 
10     CONTINUE 

COUNT=COUNT+l 

IF  (COUNT. LT. 3)  GO  TO  20 

WRITE  (5,70) 

GO  TO  50 
20     CONTINUE 

READ  (5,x,END=40,ERR=40)  IANS 

IF  (IANS)  40,30,30 
30     CONTINUE 

RETURN 
40     REWIND  5 

WRITE  (5,60) 

GO  TO  10 
50     CONTINUE 

STOP 
60     FORMAT  (IX,*  WARNING:  IMPROPER  DATA  ENTRY!   ENTER  A  POSITIVE  INTEG 

XER. ») 
70     FORMAT  (///,5X,'  PROGRAM  TERMINATION  -  2  IMPROPER  DATA  ENTRIES!') 

END 
C 
C 
C 
C 
C 
C 
c 
c 
c 
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c 

c  ==================================================================== 

C   SUBROUTINE  READC  —  INTERACTIVELY  READS  A  CHAR  STRING  REPLY. 

C   ('YES'  OR  •NO1).  IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING 

C   A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 

C  ==================================================================== 

c 
c 

SUBROUTINE  READC  (CANS) 

INTEGER  COUNT 

CHARACTERS'*  CANS 

COUNT=0 
10     CONTINUE 

COUNT=COUNT+l 

IF  (COUNT. LT. 3)  GO  TO  20 

WRITE  (5,60) 

GO  TO  40 
20     CONTINUE 

REWIND  5 

READ  (5,70,END=30,ERR=30)  CANS 

RETURN 
30     REWIND  5 

WRITE  (5,50) 

GO  TO  10 
40     CONTINUE 

STOP 
50     FORMAT  (IX,*  WARNING:   NULL  STRINGS  ARE  NOT  ALLOWED  •) 
60     FORMAT  (///,5X,«  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 
70     FORMAT  (A2) 

END 
C 
C 
C  ==================================================================== 

C   SUBROUTINE  READL  --  INTERACTIVELY  READS  A  STRING  OF  CHARACTERS. 

C   IF  THE  USER  INADVERTENTLY  ENTERS  A  NULL  STRING 

C   A  WARNING  IS  ISSUED  AND  ONE  RECOVERY  IS  ALLOWED. 

C  ==================================================================== 

c 
c 

SUBROUTINE  READL ( LLINES) 

INTEGER  COUNT,  I,  NIX 

CHARACTER**  BBLANK/1     •/,  LLINES(8) 

DO  10  1=1,8 

LLINES(I)  =  BBLANK 
10   CONTINUE 

COUNT=0 
20   COUNT=COUNT+l 

IF(C0UNT.LT.3)  GO  TO  30 

WRITE(6,70) 
GO  TO  50 
30   CONTINUE 
REWIND  5 

READ(5,80,END=40,ERR=30)(LLINES(J),J=1,9) 
RETURN 
40   REWIND  5 

WRITE(6,60) 
GO  TO  20 
50   CONTINUE 

STOP 
60   F0RMAT(1X,»  WARNING:  NULL  STRINGS  ARE  NOT  ALLOWED,  ENTER  CHARACTER 

X  VALUES. ') 
70   FORMAT  (///,5X,'  PROGRAM  TERMINATION  -  TWO  NULL  STRINGS  ENTERED!') 
80   F0RMAT(9A4) 
END 
C 
C 
C 
C 
C 
C 
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c 

c  ================================================================ 

C   SUBROUTINE  META  —  BY  CALLING  COMPRS  AS  A  SUBROUTINE  HERE, 

C   DONEPL  HAS  SUFFICIENT  TIME  TO  FINISH;  OTHERWISE  COMPRS  IGNORED 

C   AND  METAFILE  NOT  SAVED 

C  ================================================================ 

c 
c 

SUBROUTINE  META 
DO  10  1=1,900000 
10   CONTINUE 

CALL  COMPRS 

RETURN 

END 
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